Preparation

Libs and parameters

# Libraries
library(ggplot2)
# install via devtools:::install_github("JAQquent/assortedRFunctions", upgrade = "never")
library(assortedRFunctions) 
library(cowplot) 
library(stringr)
library(plyr)
library(mgcv)
library(foreach)
library(doParallel)
library(BayesFactor)
library(ciftiTools)
library(viridis)
library(tidybayes)
library(bayesplot)
library(brms)
library(data.table)
library(knitr)
library(readxl)
library(ggridges)
library(ggforce)
library(marginaleffects)
library(gghalves)
library(effectsize)
library(rstan)
library(latex2exp)

# For SVM
library(caTools) 
library(e1071) 

# Load ciftiTools and set workbench paths
possible_wb_paths <- c("/usr/bin/wb_command", "/home1/Jaquent/Toolboxes/workbench/bin_rh_linux64/")
load_ciftiTools(possible_wb_paths)
Some settings that change from computer to computer that I work on.
# Use correct locations and other settings based on computer
if(Sys.info()[4] == "DESKTOP-335I26I"){
  # Work laptop (Windows)
  path2imaging_results2 <- "D:/Seafile/imaging_results"
} else if(Sys.info()[4] == 'DESKTOP-91CQCSQ') {
  # Work desktop (Windows)
  path2imaging_results2 <- "D:/imaging_results"
} else if(Sys.info()[4] == 'alex-Zenbook-UX3404VA-UX3404VA') {
  # Work laptop (Linux)
  path2imaging_results2 <- "/media/alex/shared/Seafile/imaging_results"
} else if(Sys.info()[4] == "greengoblin"){
  # Main desktop PC (Linux)
  path2imaging_results2 <- "/media/alex/work/Seafile/imaging_results" 
} else if(Sys.info()[4] == "GREEN-GOBLIN-WI"){
  # Main desktop PC 
  path2imaging_results2 <- "E:/Seafile/imaging_results" 
} else {
  # Personal laptop (Windows)
  path2imaging_results2 <- "D:/OLM/imaging_results"
}

# Seed
set.seed(20240205)
Click here for detailed session information.
sessioninfo::session_info()
## ─ Session info ───────────────────────────────────────────────────────────────
##  setting  value
##  version  R version 4.4.0 (2024-04-24)
##  os       Ubuntu 22.04.5 LTS
##  system   x86_64, linux-gnu
##  ui       X11
##  language en_GB:en
##  collate  en_GB.UTF-8
##  ctype    en_GB.UTF-8
##  tz       Asia/Shanghai
##  date     2025-10-30
##  pandoc   3.1.1 @ /usr/lib/rstudio/resources/app/bin/quarto/bin/tools/ (via rmarkdown)
## 
## ─ Packages ───────────────────────────────────────────────────────────────────
##  package            * version    date (UTC) lib source
##  abind                1.4-5      2016-07-21 [1] CRAN (R 4.4.0)
##  arrayhelpers         1.1-0      2020-02-04 [1] CRAN (R 4.4.0)
##  assortedRFunctions * 1.1.3      2025-09-28 [1] Github (JAquent/assortedRFunctions@a84af63)
##  backports            1.5.0      2024-05-23 [1] CRAN (R 4.4.0)
##  base64enc            0.1-3      2015-07-28 [1] CRAN (R 4.4.0)
##  BayesFactor        * 0.9.12-4.7 2024-01-24 [1] CRAN (R 4.4.0)
##  bayesplot          * 1.11.1     2024-02-15 [1] CRAN (R 4.4.0)
##  bayestestR           0.16.1     2025-07-01 [1] CRAN (R 4.4.0)
##  bitops               1.0-7      2021-04-24 [1] CRAN (R 4.4.0)
##  bridgesampling       1.1-2      2021-04-16 [1] CRAN (R 4.4.0)
##  brms               * 2.22.0     2024-09-23 [1] CRAN (R 4.4.0)
##  Brobdingnag          1.2-9      2022-10-19 [1] CRAN (R 4.4.0)
##  bslib                0.7.0      2024-03-29 [1] CRAN (R 4.4.0)
##  cachem               1.1.0      2024-05-16 [1] CRAN (R 4.4.0)
##  caTools            * 1.18.2     2021-03-28 [1] CRAN (R 4.4.0)
##  cellranger           1.1.0      2016-07-27 [1] CRAN (R 4.4.0)
##  checkmate            2.3.1      2023-12-04 [1] CRAN (R 4.4.0)
##  ciftiTools         * 0.17.4     2025-01-23 [1] CRAN (R 4.4.0)
##  class                7.3-20     2022-01-13 [4] CRAN (R 4.1.2)
##  cli                  3.6.2      2023-12-11 [1] CRAN (R 4.4.0)
##  coda               * 0.19-4.1   2024-01-31 [1] CRAN (R 4.4.0)
##  codetools            0.2-18     2020-11-04 [4] CRAN (R 4.0.3)
##  colorspace           2.1-0      2023-01-23 [1] CRAN (R 4.4.0)
##  cowplot            * 1.1.3      2024-01-22 [1] CRAN (R 4.4.0)
##  data.table         * 1.17.8     2025-07-10 [1] CRAN (R 4.4.0)
##  datawizard           1.1.0      2025-05-09 [1] CRAN (R 4.4.0)
##  digest               0.6.35     2024-03-11 [1] CRAN (R 4.4.0)
##  distributional       0.4.0      2024-02-07 [1] CRAN (R 4.4.0)
##  doParallel         * 1.0.17     2022-02-07 [1] CRAN (R 4.4.0)
##  dplyr                1.1.4      2023-11-17 [1] CRAN (R 4.4.0)
##  e1071              * 1.7-14     2023-12-06 [1] CRAN (R 4.4.0)
##  effectsize         * 0.8.8      2024-05-12 [1] CRAN (R 4.4.0)
##  emmeans              1.10.4     2024-08-21 [1] CRAN (R 4.4.0)
##  estimability         1.5.1      2024-05-12 [1] CRAN (R 4.4.0)
##  evaluate             0.23       2023-11-01 [1] CRAN (R 4.4.0)
##  fansi                1.0.6      2023-12-08 [1] CRAN (R 4.4.0)
##  farver               2.1.2      2024-05-13 [1] CRAN (R 4.4.0)
##  fastmap              1.2.0      2024-05-15 [1] CRAN (R 4.4.0)
##  foreach            * 1.5.2      2022-02-02 [1] CRAN (R 4.4.0)
##  generics             0.1.3      2022-07-05 [1] CRAN (R 4.4.0)
##  ggdist               3.3.2      2024-03-05 [1] CRAN (R 4.4.0)
##  ggforce            * 0.4.2      2024-02-19 [1] CRAN (R 4.4.0)
##  gghalves           * 0.1.4      2022-11-20 [1] CRAN (R 4.4.0)
##  ggplot2            * 3.5.2      2025-04-09 [1] CRAN (R 4.4.0)
##  ggridges           * 0.5.6      2024-01-23 [1] CRAN (R 4.4.0)
##  gifti                0.8.0      2020-11-11 [1] CRAN (R 4.4.0)
##  glue                 1.7.0      2024-01-09 [1] CRAN (R 4.4.0)
##  gridExtra            2.3        2017-09-09 [1] CRAN (R 4.4.0)
##  gtable               0.3.5      2024-04-22 [1] CRAN (R 4.4.0)
##  htmltools            0.5.8.1    2024-04-04 [1] CRAN (R 4.4.0)
##  inline               0.3.19     2021-05-31 [1] CRAN (R 4.4.0)
##  insight              1.3.1      2025-06-30 [1] CRAN (R 4.4.0)
##  iterators          * 1.0.14     2022-02-05 [1] CRAN (R 4.4.0)
##  jquerylib            0.1.4      2021-04-26 [1] CRAN (R 4.4.0)
##  jsonlite             1.8.8      2023-12-04 [1] CRAN (R 4.4.0)
##  knitr              * 1.50       2025-03-16 [1] CRAN (R 4.4.0)
##  latex2exp          * 0.9.6      2022-11-28 [1] CRAN (R 4.4.0)
##  lattice              0.20-45    2021-09-22 [4] CRAN (R 4.1.1)
##  lifecycle            1.0.4      2023-11-07 [1] CRAN (R 4.4.0)
##  loo                  2.8.0      2024-07-03 [1] CRAN (R 4.4.0)
##  magrittr             2.0.3      2022-03-30 [1] CRAN (R 4.4.0)
##  marginaleffects    * 0.28.0     2025-06-25 [1] CRAN (R 4.4.0)
##  MASS                 7.3-64     2025-01-04 [1] CRAN (R 4.4.0)
##  Matrix             * 1.7-0      2024-04-26 [1] CRAN (R 4.4.0)
##  MatrixModels         0.5-3      2023-11-06 [1] CRAN (R 4.4.0)
##  matrixStats          1.3.0      2024-04-11 [1] CRAN (R 4.4.0)
##  mgcv               * 1.8-39     2022-02-24 [4] CRAN (R 4.1.2)
##  multcomp             1.4-26     2024-07-18 [1] CRAN (R 4.4.0)
##  munsell              0.5.1      2024-04-01 [1] CRAN (R 4.4.0)
##  mvtnorm              1.2-5      2024-05-21 [1] CRAN (R 4.4.0)
##  nlme               * 3.1-155    2022-01-13 [4] CRAN (R 4.1.2)
##  oro.nifti            0.11.4     2022-08-10 [1] CRAN (R 4.4.0)
##  parameters           0.27.0     2025-07-09 [1] CRAN (R 4.4.0)
##  pbapply              1.7-2      2023-06-27 [1] CRAN (R 4.4.0)
##  pillar               1.9.0      2023-03-22 [1] CRAN (R 4.4.0)
##  pkgbuild             1.4.4      2024-03-17 [1] CRAN (R 4.4.0)
##  pkgconfig            2.0.3      2019-09-22 [1] CRAN (R 4.4.0)
##  plyr               * 1.8.9      2023-10-02 [1] CRAN (R 4.4.0)
##  polyclip             1.10-6     2023-09-27 [1] CRAN (R 4.4.0)
##  posterior            1.6.1      2025-02-27 [1] CRAN (R 4.4.0)
##  proxy                0.4-27     2022-06-09 [1] CRAN (R 4.4.0)
##  purrr                1.0.2      2023-08-10 [1] CRAN (R 4.4.0)
##  QuickJSR             1.2.2      2024-06-07 [1] CRAN (R 4.4.0)
##  R.methodsS3          1.8.2      2022-06-13 [1] CRAN (R 4.4.0)
##  R.oo                 1.26.0     2024-01-24 [1] CRAN (R 4.4.0)
##  R.utils              2.12.3     2023-11-18 [1] CRAN (R 4.4.0)
##  R6                   2.5.1      2021-08-19 [1] CRAN (R 4.4.0)
##  RColorBrewer         1.1-3      2022-04-03 [1] CRAN (R 4.4.0)
##  Rcpp               * 1.0.12     2024-01-09 [1] CRAN (R 4.4.0)
##  RcppParallel         5.1.7      2023-02-27 [1] CRAN (R 4.4.0)
##  readxl             * 1.4.3      2023-07-06 [1] CRAN (R 4.4.0)
##  rlang                1.1.4      2024-06-04 [1] CRAN (R 4.4.0)
##  rmarkdown            2.27       2024-05-17 [1] CRAN (R 4.4.0)
##  RNifti               1.6.1      2024-03-07 [1] CRAN (R 4.4.0)
##  rstan              * 2.32.6     2024-03-05 [1] CRAN (R 4.4.0)
##  rstantools           2.4.0      2024-01-31 [1] CRAN (R 4.4.0)
##  rstudioapi           0.16.0     2024-03-24 [1] CRAN (R 4.4.0)
##  sandwich             3.1-1      2024-09-15 [1] CRAN (R 4.4.0)
##  sass                 0.4.9      2024-03-15 [1] CRAN (R 4.4.0)
##  scales               1.3.0      2023-11-28 [1] CRAN (R 4.4.0)
##  sessioninfo          1.2.2      2021-12-06 [1] CRAN (R 4.4.0)
##  StanHeaders        * 2.32.9     2024-05-29 [1] CRAN (R 4.4.0)
##  stringi              1.8.4      2024-05-06 [1] CRAN (R 4.4.0)
##  stringr            * 1.5.1      2023-11-14 [1] CRAN (R 4.4.0)
##  survival             3.2-13     2021-08-24 [4] CRAN (R 4.1.1)
##  svUnit               1.0.6      2021-04-19 [1] CRAN (R 4.4.0)
##  tensorA              0.36.2.1   2023-12-13 [1] CRAN (R 4.4.0)
##  TH.data              1.1-2      2023-04-17 [1] CRAN (R 4.4.0)
##  tibble               3.2.1      2023-03-20 [1] CRAN (R 4.4.0)
##  tidybayes          * 3.0.6      2023-08-12 [1] CRAN (R 4.4.0)
##  tidyr                1.3.1      2024-01-24 [1] CRAN (R 4.4.0)
##  tidyselect           1.2.1      2024-03-11 [1] CRAN (R 4.4.0)
##  tweenr               2.0.3      2024-02-26 [1] CRAN (R 4.4.0)
##  utf8                 1.2.4      2023-10-22 [1] CRAN (R 4.4.0)
##  vctrs                0.6.5      2023-12-01 [1] CRAN (R 4.4.0)
##  viridis            * 0.6.5      2024-01-29 [1] CRAN (R 4.4.0)
##  viridisLite        * 0.4.2      2023-05-02 [1] CRAN (R 4.4.0)
##  withr                3.0.0      2024-01-16 [1] CRAN (R 4.4.0)
##  xfun                 0.52       2025-04-02 [1] CRAN (R 4.4.0)
##  xml2                 1.3.6      2023-12-04 [1] CRAN (R 4.4.0)
##  xtable               1.8-4      2019-04-21 [1] CRAN (R 4.4.0)
##  yaml                 2.3.8      2023-12-11 [1] CRAN (R 4.4.0)
##  zoo                  1.8-12     2023-04-13 [1] CRAN (R 4.4.0)
## 
##  [1] /home/alex/R/x86_64-pc-linux-gnu-library/4.4
##  [2] /usr/local/lib/R/site-library
##  [3] /usr/lib/R/site-library
##  [4] /usr/lib/R/library
## 
## ──────────────────────────────────────────────────────────────────────────────
Click here for chunk for figure and other reporting parameters.
# Significance cut off 
cutOff <- 1.301 # equals -log(0.05, 10)

# Parameters how to report means
report_type   <- 1
digits1       <- 2
rounding_type <- "signif"

theme_journal <- function(base_size = 7, linewidth = 0.35, theme_name = "classic") {
  # Select a base theme from https://ggplot2.tidyverse.org/reference/ggtheme.html
  if(theme_name == "grey" | theme_name == "gray"){
    base_theme <- theme_grey(base_size = base_size)
    
  } else if(theme_name == "bw"){
    base_theme <- theme_bw(base_size = base_size)
    
  } else if(theme_name == "linedraw"){
    base_theme <- theme_linedraw(base_size = base_size)
    
  } else if(theme_name == "light"){
    base_theme <- theme_light(base_size = base_size)
    
  } else if(theme_name == "dark"){
    base_theme <- theme_dark(base_size = base_size)
    
  } else if(theme_name == "minimal"){
    base_theme <- theme_minimal(base_size = base_size)
    
  } else if(theme_name == "classic"){
    base_theme <- theme_classic(base_size = base_size)
    
  } else {
    stop("Unknown theme. Check https://ggplot2.tidyverse.org/reference/ggtheme.html")
  }
  
  # Base them plus customisation
  base_theme +
    theme(text = element_text(color = "black", size = base_size),
          axis.text = element_text(size = base_size * 0.9),  # 90% of base
          axis.line = element_line(linewidth = linewidth),
          axis.ticks = element_line(linewidth = linewidth),
          legend.text = element_text(size = base_size * 0.9), # 90% of base
          legend.key.size = unit(0.5, "lines"),
          plot.title = element_text(hjust = 0.5, size = base_size))
}

# Updating geom's defaults
update_geom_defaults("point", list(size = 0.5, stroke = 0.25))
update_geom_defaults("line", list(linewidth = 0.10))
update_geom_defaults("density", list(linewidth = 0.25))
update_geom_defaults("boxplot", list(linewidth = 0.25, size = 0.2))
update_geom_defaults("smooth", list(linewidth = 0.5))
update_geom_defaults("hline", list(linewidth = 0.25))
update_geom_defaults("vline", list(linewidth = 0.25))
update_geom_defaults("curve", list(linewidth = 0.25))
update_geom_defaults("segment", list(linewidth = 0.25))

# Parameters for saving figures
## Information: https://www.nature.com/ncomms/submit/how-to-submit
### PDF page: 210 x 276 mm
### Column widths in mm
single_column <- 88
double_column <- 180
dpi           <- 1000
figurePath    <- "figures/SpaNov/"

# Colours used for visualisation
boxplot_border_col  <- "black"
boxplot_point_col   <- "darkgrey"
boxplot_line_col    <- "darkgrey"
base_col            <- c("#7091C0", "#4A66AC", "#364875")
spectral_colours    <- c("#5956a5", "#a60a44")
cool_warm_colours   <- c("#3C4DC1", "#B70B28")
novFam_gradient     <- viridis(n = 6, option = "H", direction = -1)

# Information for Yeo 7
## Names etc. for the networks in Yeo 7
Yeo7_fullNames <- c("Frontoparietal", "Default", "Dorsal Attention", "Limbic", 
                    "Ventral Attention", "Somatomotor", "Visual")
Yeo7_abbr      <- c("Cont", "Default", "DorsAttn", "Limbic", "SalVentAttn", "SomMot", "Vis")

## Colours used for Yeo 7
Yeo7_colours <- c("#E79523", "#CD3E4E", "#00760F", "#DCF8A4", "#C43BFA", "#4682B4", "#781286")

Cifti files

Click here for support CIFTI files etc.
# Loading the support CIFTI files
## Place where to find some of the CIFTI files and parcellations
CIFTI_locations <- "data/ignore_fMRI_version1/sourceFiles/"

## CIFTI files
parcellationFile <- "Q1-Q6_RelatedValidation210.CorticalAreas_dil_Final_Final_Areas_Group_Colors_with_Atlas_ROIs2.32k_fs_LR.dlabel.nii"
CAB_NP           <- "CortexSubcortex_ColeAnticevic_NetPartition_wSubcorGSR_netassignments_LR.dlabel.nii"
surfLeft         <- "S1200.L.inflated_MSMAll.32k_fs_LR.surf.gii"
surfRight        <- "S1200.R.inflated_MSMAll.32k_fs_LR.surf.gii"

## Combine CIFTI_locations with file name
parcellationFile <- paste0(CIFTI_locations, parcellationFile)
CAB_NP           <- paste0(CIFTI_locations, CAB_NP)
surfLeft         <- paste0(CIFTI_locations, surfLeft)
surfRight        <- paste0(CIFTI_locations, surfRight)

## Loading CIFTIw via ciftiTools as xiis
### Get MMP parcellation
MMP_xii <- ciftiTools::read_cifti(parcellationFile,
                                  brainstructures = "all", 
                                  surfL_fname = surfLeft, 
                                  surfR_fname = surfRight)

### Load Yeo 7 parcellation
Yeo7_xii <- ciftiTools::read_cifti("other_stuff/Yeo7.dlabel.nii",
                                  surfL_fname = surfLeft, 
                                  surfR_fname = surfRight)

## Load other stuff
### Load the parcel names for MMP
parcel_names <- read.csv("data/ignore_fMRI_version1/extracted_values/Parcellations/MP1.0_210V_parcel_names.csv", header = FALSE)

### Load the extracted MNI coordinates
MNI_coord <- read.csv("other_stuff/cifti_subcortical_MNI152_coordinates.csv")

### Load the hippocampal projection values
load("other_stuff/projected_HC.RData")

Custom functions

Click here for code for custom functions.
create_str_from_avg_comparisons <- function (x, measure_name = "X" , digits = 2, rounding_type = "round"){
  # Convert to numeric values
  x_num <- suppressWarnings(as.numeric(x))
  
  # Round and create string
  if (rounding_type == "signif") {
      return_string <- paste0(signif(x_num[3], digits), " ", measure_name ,  " (95 % CI [", 
          signif(x_num[4], digits), ", ", signif(x_num[5], digits), 
          "])")
  }
  else if (rounding_type == "round") {
      return_string <- paste0(round(x_num[3], digits), " ", measure_name ,  " (95 % CI [", 
          round(x_num[4], digits), ", ", round(x_num[5], digits), 
          "])")
  }
  else {
      stop("Wrong rounding type. Choose signif or round.")
  }
  return(return_string)
}

Loading and preparing the data

Click here for loading data.
# Load the data
## Specify paths where the data is saved
path2data <- "data/ignore_fMRI_version1/"
EV_folder <- "ignore_eventTable2/"

## Load the look-up table that contains information of R-numbers which are retracted 
lookupTable  <- read.csv(paste0(path2data, "lookUpTable.csv"))

## Load .RData images of the combined data (all subject in one DF)
load(paste0(path2data, "combined_data/demographics.RData"))
load(paste0(path2data, "combined_data/DW_all_data.RData"))
load(paste0(path2data, "combined_data/OLM_7T_all_data.RData"))
load(paste0(path2data, "combined_data/OLM_3T_all_data.RData"))
load(paste0(path2data, "combined_data/question_data.RData"))

# Select the subjects included in this analysis
## Load the subjects that are included in this analysis
subjectFile <- readLines(paste0(path2data, "SpaNov_subject2analyse.txt"))
subjIDs_R   <- str_split(subjectFile, pattern = ",")[[1]] 
subjIDs     <- lookupTable$anonKey[lookupTable$Rnum %in% subjIDs_R]
# Important note: subjIDs_R do not have the same order as subjIDs!!!!!!!!!!!!!

## Subset to data that is being included in the analysis
OLM_7T_position_data <- OLM_7T_position_data[OLM_7T_position_data$subject %in% subjIDs, ]
demographics         <- demographics[demographics$subject %in% subjIDs, ]
DW_position_data     <- DW_position_data[DW_position_data$subject %in% subjIDs, ]
OLM_7T_logEntries    <- OLM_7T_logEntries[OLM_7T_logEntries$ppid %in% subjIDs, ]
OLM_7T_trial_results <- OLM_7T_trial_results[OLM_7T_trial_results$subject %in% subjIDs, ]
question_data        <- question_data[question_data$subject %in% subjIDs, ]

# Create positionData
positionData <- OLM_7T_position_data

# Add ppid to this data frame
subjects_anon     <- unique(positionData$subject)
positionData$ppid <- NA
for(i in 1:length(subjects_anon)){
  positionData$ppid[positionData$subject == subjects_anon[i]] <- lookupTable$Rnum[lookupTable$anonKey == subjects_anon[i]]
}

# Add ppid to OLM_7T_trial_results
OLM_7T_trial_results$ppid <- NA
for(i in 1:nrow(OLM_7T_trial_results)){
  OLM_7T_trial_results$ppid[i] <- lookupTable$Rnum[lookupTable$anonKey == OLM_7T_trial_results$subject[i]]
}

# Subset to retrieval only
OLM_7T_retrieval <- OLM_7T_trial_results[OLM_7T_trial_results$trialType == "retrieval", ]

# Get the object locations to verify object placement in screenshots
obj_locations <- ddply(OLM_7T_trial_results, c("targets", "objectName", "targetNames"),
                       summarise, object_x_sd = sd(object_x), object_x = mean(object_x),
                       object_z_sd = sd(object_z), object_z = mean(object_z))

# Subset to encoding only
OLM_7T_encoding <- OLM_7T_trial_results[OLM_7T_trial_results$trialType == "encoding", ]

# Get cue times
cues_agg <- data.frame(ppid = OLM_7T_encoding$ppid,
                       trial = OLM_7T_encoding$trial_num,
                       start = OLM_7T_encoding$start_time,
                       end = OLM_7T_encoding$start_time + OLM_7T_encoding$cue,
                       duration = OLM_7T_encoding$cue,
                       block_num = OLM_7T_encoding$block_num)
# Get delay times
delays_agg <- data.frame(ppid = OLM_7T_encoding$ppid,
                         trial = OLM_7T_encoding$trial_num,
                         start = OLM_7T_encoding$start_time + OLM_7T_encoding$cue,
                         end = OLM_7T_encoding$start_time + OLM_7T_encoding$cue + OLM_7T_encoding$delay,
                         duration = OLM_7T_encoding$delay,
                         block_num = OLM_7T_encoding$block_num)
# Parameter
downsample_value  <- 0.2 #200 msec/0.2 sec

# Add run information
run_info <- ddply(OLM_7T_trial_results, c("block_num", "trial_num"), summarise, N = length(ppid))
positionData$run <- NA

# Loop through all trials
for(i in 1:nrow(run_info)){
  positionData$run[positionData$trial == run_info$trial_num[i]] <- run_info$block_num[i]
}

# Prepare the cluster for a parallel loop
# Create the cluster following https://www.blasbenito.com/post/02_parallelizing_loops_with_r/
my.cluster <- parallel::makeCluster(detectCores() - 2, type = "PSOCK")

# Register it to be used by %dopar%
doParallel::registerDoParallel(cl = my.cluster)

# Algorithm to down sample to a sample every x msec
# For that loop through tempData_currentSubj and check if 200 msec have passed since current time
# Each row where a new run begins is always included. After that only a sample that is 200 msec
# passed the current time. Especially the first sample of a run is important to accurately
# determine when the first image is recorded and that all the onsets correspond to that. 

# Extract unique participants but in order in which they appear. This was verified
# E.g. via positionData$ppid[!duplicated(positionData$ppid)]
subjects <- unique(positionData$ppid)

# Run parallel loop
include <- foreach(i = 1:length(subjIDs_R), .combine = 'c') %dopar% {
  # Subset to current subject
  tempData_currentSubj <- positionData[positionData$ppid == subjIDs_R[i], ]
  
  # Create include variable
  tempInclude <- rep(FALSE, nrow(tempData_currentSubj))
  
  # Loop through the each row
  for(i in 1:nrow(tempData_currentSubj)){
    # First time
    if(i == 1){
      currentTime <- tempData_currentSubj$time[i]
      tempInclude[i]  <- TRUE
    } else {
      # Check if the last row was from a different run
      if(tempData_currentSubj$run[i - 1] != tempData_currentSubj$run[i]){
        # Set new time in this case
        currentTime <- tempData_currentSubj$time[i]
        tempInclude[i]  <- TRUE
      } else {
        # Check ih ith time is larger than currentTime + downsample_value
        if(tempData_currentSubj$time[i] > currentTime + downsample_value){
          # Set new time in this case
          currentTime <- tempData_currentSubj$time[i]
          tempInclude[i] <- TRUE
        }
      }
    }
  }
  # Return
  tempInclude
}

# Stop cluster again
parallel::stopCluster(cl = my.cluster)

# Apply downsampling to positionData
positionData2 <- positionData[include,]


# Check if there is no problem
time_diff <- ddply(positionData2, c("ppid", "run", "trial"), summarise, 
                   mean_time_diff = mean(diff(time)),
                   sd_time_diff   = sd(diff(time)))
# Now I manually checked the run start time of the first two participants

# Get the number of seeds of this analysis
numSeeds  <- 10
limValues <- c(-90,90)

# Bin the environment
x           <- positionData2$pos_x
z           <- positionData2$pos_z
env_sectors <- voronoi_tessellation_grid_binning_2d(x, z,limValues, numSeeds, "hexagon", 
                                                useParallelisation  = TRUE)

# Add result back to the df (sector2 is only for plotting)
positionData2$sector  <- env_sectors
positionData2$sector2 <- factor(env_sectors, levels = sample(unique(env_sectors)))

# Save in intermediate data
save(positionData2, file = "intermediate_data/positionData2.RData")
# Load data
load("intermediate_data/positionData2.RData")

# Exclude control trials
no_control <- positionData2[positionData2$trialType != "control", ]

# Include only Run 1, 2 and 3 because the last retrieval run 
# doesn't count as far as this analysis is concerned
no_control <- no_control[no_control$run %in% 1:3, ]

Results

Modelling of spatial novelty during naturalistic navigation

Reporting demographic information

n       <- nrow(demographics)
str1    <- mean_SD_str2(demographics$age, type = 1, digits = digits1, rounding_type = rounding_type, measure = "years")
females <- table(demographics$gender)[1]
males   <- table(demographics$gender)[2]
  • n = 56
  • Age = M = 24 years (SD = 2.8 years)
  • Gender ratio = 35/21
  • Age range 19.7, 37.3

Figure 1

3D histogramms

get_seeds_from_voronoi_tessellation_hexagon <- function(limValues, numSeeds){
  # Create bins like it is done in voronoi_tessellation_grid_binning_2d
  # This is code directly from that function to get the coordinates of each sector
  xLim    <- sort(limValues)
  yLim    <- sort(limValues)
  byValue <- ((xLim[2] - xLim[1])/(numSeeds - 1))
  xRange  <- seq(from = xLim[1], to = xLim[2], by = byValue)
  xStepSize <- xRange[2] - xRange[1]
  x_shiftValue <- xStepSize/4
  yRange <- seq(from = yLim[1], to = yLim[2], by = xStepSize * (sqrt(3)/2))
  yRange <- yRange[1:length(xRange)]
  yRange <- yRange + (yLim[2] - max(yRange))/2
  seeds <- data.frame(x = xRange[1], y = yRange)
  shiftIndex <- rep(c(1, -1), length.out = nrow(seeds))
  seeds$x <- seeds$x + shiftIndex * x_shiftValue
  for(i in 2:length(xRange)){
    seeds <- rbind(seeds, data.frame(x = xRange[i] + shiftIndex * x_shiftValue, y = yRange))
  }
  seeds$sector <- row.names(seeds)
  return(seeds)
}

This creates the data necessary to create the 3D histogram using blender.

# Load data from 
load("E:/research_projects/OLM_project/analysis/ignore_eventTable3/images/SpaNov_event_file_contPM.RData")

# Convert list to data frame
data <- rbindlist(subj_list, idcol = "subject")

# Subject to character
data$subject <- as.character(data$subject)

# Subset to run 1, 2, 3 to include the first retrieval run
data_sub2 <- data[data$run != 4, ]

# Calculate average/sum per subject
sector_average <- ddply(data_sub2, c("subject","sector"), summarise, 
                        visits = sum(visits, na.rm = TRUE),
                        lastVisit = mean(lastVisit, na.rm = TRUE))

# Calculate average per sector
sector_average <- ddply(sector_average, c("sector"), summarise, 
                        visits = mean(visits, na.rm = TRUE),
                        lastVisit = mean(lastVisit, na.rm = TRUE))

# Bin the environment
limValues   <- c(-90,90)
numSeeds    <- 10 # Number of values per axis

# Create bins like it is done in voronoi_tessellation_grid_binning_2d
seeds <- get_seeds_from_voronoi_tessellation_hexagon(limValues, numSeeds)

# Add x & y coordinates to sector_visits. For this just loop the df
sector_average$x <- NA
sector_average$y <- NA
for(i in 1:nrow(sector_average)){
  # Get current sector
  currentSector <- sector_average$sector[i]
  
  # Get corresponding row index in seeds
  rowID <- which(seeds$sector == currentSector)
  
  # Assign the correct coordinates based on rowID
  sector_average$x[i] <- seeds$x[rowID]
  sector_average$y[i] <- seeds$y[rowID]
}


# Prepare data for 3D plot for the average number of visits per persons to each sector
temp_data <- sector_average[,-3]

# Write .csv file
write.csv(x = temp_data, file = "other_stuff/source_data_files/blenderData_avg_visits.csv", 
          quote = FALSE, row.names = FALSE)

# Prepare data for 3D plot for the average number of visits per persons to each sector
## Remove NA values & remove wrong column
temp_data <- na.omit(sector_average)
temp_data <- temp_data[,-2]

# Write .csv file
write.csv(x = temp_data, file = "other_stuff/source_data_files/blenderData_avg_time_since_last_visit.csv", 
          quote = FALSE, row.names = FALSE)

# Colours
barColours <- viridisLite::viridis(n = 5, option = "D")

Range of values is 75 to 111321.

Distribution of total path lengths

# Subset to only data from the encoding part
pathData <- positionData[positionData$trialType == 'encoding', ]

# Calculate the how much the perfect participant would need to travel
## Subset to only encoding trials
boolIndex       <- OLM_7T_trial_results$trialType == "encoding"
OLM_7T_encoding <- OLM_7T_trial_results[boolIndex, ]

# Function to calculate the path lengths 
calculate_path_length <- function(pos_x, pos_z){
  # Calculate difference between points and square them
  diff_X <- diff(pos_x)^2
  diff_z <- diff(pos_z)^2
  
  # Calculate the sum and take the square root
  dists_travelled <- sqrt(diff_X + diff_z)
  
  # Sum up all distances travelled to one path_length
  return(sum(dists_travelled))
}

# Calculate the path traveled for each trial and each subject
# Use pathData because it only includes encoding data
path_lenghts <- ddply(pathData, c("ppid", "trial"), 
                      summarise, 
                      distance = calculate_path_length(pos_x, pos_z))

# Calculate the total path lengths for the whole experiment
total_path_lengths <- ddply(path_lenghts, c("ppid"), 
                            summarise, distance = sum(distance))

# Write .csv file
write.csv(x = total_path_lengths, file = "other_stuff/source_data_files/total_path_lengths.csv", 
          quote = FALSE, row.names = FALSE)

# Create a box plot 
total_path_lengths_plot <- ggplot(total_path_lengths, aes(x = distance)) + 
  geom_density(fill = base_col[1]) + 
  geom_jitter(aes(y = -0.001), height = 0.0005, width = 0, colour = boxplot_point_col) +
  labs(y = "Density", x = "Total path length (vm)", title = "") +
  scale_x_continuous(breaks = c(3000, 3300, 3600)) + 
  coord_cartesian(xlim = c(3000, 3600), expand = FALSE, ylim = c(-0.001*2, 0.0065)) +
  theme_journal() +
  theme(axis.text.y = element_blank(), 
        axis.ticks.y = element_blank(),
        axis.line.y = element_blank(),
        plot.margin = unit(c(1, 2.7, 1, 2.7), "mm"))

Distributions of event durations

# Load the event files
event_files <- list.files(path = "event_tables/OLMe_7T_SpaNov_gradient_6lvl/", recursive = TRUE, full.names = TRUE)

# Remove "per_tra.txt" and "dur_tra.txt" from list
event_files <- event_files[!str_detect(event_files, "tra.txt")]

# Remove cue and delay
event_files <- event_files[!str_detect(event_files, "cue.txt") & !str_detect(event_files, "delay.txt")]

# Create list
event_list <- list()

# Load the event files
for(i in 1:length(event_files)){
  # Load current event file
  temp_event_file        <- read.table(event_files[i], header = FALSE, sep = "\t")
  temp_event_file$V3     <- NULL # Remove the third column
  names(temp_event_file) <- c("onset", "duration")
  
  # Split the file path to get information
  event_file_path_split   <- str_split_1(event_files[i], pattern = "/")
  temp_event_file$subject <- event_file_path_split[4]
  temp_event_file$run     <- event_file_path_split[5]
  
  # To extract the level of the event extract "lvl" plus an integer from string
  temp_event_file$lvl     <- str_extract(event_file_path_split[8], "lvl[0-9]+")
  
  # Add to list
  event_list[[i]] <- temp_event_file
}

# Combine to data frame
event_df <- as.data.frame(rbindlist(event_list))

# Calculate offset
event_df$offset <- event_df$onset + event_df$duration

# Convert lvl label to number
event_df$lvl_num <- as.numeric(str_extract(event_df$lvl, "[0-9]+"))

# Change run label
event_df$run <- ifelse(event_df$run == "run-01", "Run 1", "Run 2")

# Write .csv file
write.csv(x = event_df, file = "other_stuff/source_data_files/event_duration.csv", 
          quote = FALSE, row.names = FALSE)

# Create a box plot 
event_duration_plot <- ggplot(event_df, aes(x = duration)) + 
  geom_density(fill = base_col[1]) + 
  geom_jitter(aes(y = -0.1), height = 0.05, width = 0, colour = boxplot_point_col, alpha = 0.1) +
  labs(y = "Density", x = "Event duration (s)", title = "") +
  scale_x_continuous(breaks = c(0, 16, 32)) + 
  coord_cartesian(xlim = c(0, 32), expand = FALSE, ylim = c(-0.1*2, 0.8)) +
  theme_journal() +
  theme(axis.text.y = element_blank(), 
        axis.ticks.y = element_blank(),
        axis.line.y = element_blank(),
        plot.margin = unit(c(1, 2.5, 1, 2.5), "mm"))

Distributions of Rayleigh vector length

Calculate Rayleigh Test of Uniformity for each sector

#install.packages("circular")  # If not already installed
library(circular)

# Function to calculate the test and apply it to the data frame
calc_rayleigh_test_FUN <- function(angles){
  # Create circular object
  circ_data <- circular(angles, units = "degrees", template = "none")
 
 # Calculate Rayleigh test which includes R
 rayleigh_test <- rayleigh.test(circ_data)
 
 # Return
 return(c(rayleigh_test$statistic, rayleigh_test$p.value))
}

# Apply test to all sectors
sector_rayleigh <- ddply(no_control, c("sector"), summarise,
                         resultant_length = calc_rayleigh_test_FUN(rot_y)[1],
                         p.value = calc_rayleigh_test_FUN(rot_y)[2])

# Write .csv file
write.csv(x = sector_rayleigh, file = "other_stuff/source_data_files/sector_rayleigh.csv", 
          quote = FALSE, row.names = FALSE)

sector_rayleigh_plot <- ggplot(sector_rayleigh, aes(x = resultant_length)) +
  geom_histogram(colour = base_col[1], fill = base_col[1]) + 
  geom_vline(xintercept = mean(sector_rayleigh$resultant_length), colour = "black", linetype = 2, size = 0.5) +
  scale_x_continuous(breaks = c(0, 0.5, 1)) + 
  scale_y_continuous(breaks = c(0, 5, 10)) + 
  coord_cartesian(xlim = c(-0.02, 1), ylim = c(0, 10), 
                  expand = FALSE) +
  labs(x = "Length of Rayleigh vector", y = "Count") +
  theme_journal() +
  theme(plot.margin = unit(c(1, 2.5, 1, 2.5), "mm"))

# Report mean
str1 <- mean_SD_str2(sector_rayleigh$resultant_length, report_type, digits1, rounding_type)

Time spend in each locomotion state

Calculate overall time spent on translation, rotation and being stationary:

# Remove cue & delay periods
locomotion_data         <- pathData
locomotion_data$include <- TRUE
subjs  <- unique(locomotion_data$ppid)
trials <- unique(locomotion_data$trial)
for(i in 1:length(subjs)){
  for(j in 1:length(trials)){
    # Get cue start & delay end time
    cue_start <- cues_agg$start[cues_agg$ppid == subjs[i] & cues_agg$trial == trials[j]]
    delay_end <- delays_agg$end[delays_agg$ppid == subjs[i] & delays_agg$trial == trials[j]]
    
    # Set include to false for those times
    bool_index <- locomotion_data$ppid == subjs[i] & 
                  locomotion_data$trial == trials[j] &
                  locomotion_data$time >= cue_start & 
                  locomotion_data$time <= delay_end
    locomotion_data$include[bool_index] <- FALSE
  }
}

# Remove those times
locomotion_data <- locomotion_data[locomotion_data$include, ]

# The amount we found the values to avoid false positive
rotation_round  <- 2 # round rotation values to this decimal point

# Function to determine which state a time point belongs to
what_state <- function(rot_y, moving2){
  # Get angles 
  angle1 <- rot_y[2:length(rot_y)]
  angle2 <- rot_y[1:(length(rot_y)-1)]
  
  # Calculate the amount was rotated between the time points and then rotate
  rotated <- c(NA, round(angularDifference(angle1, angle2), rotation_round))
  
  # If rotation is zero called it stationary, otherwise rotation
  tra_rot_sta <- ifelse(abs(rotated) == 0 | is.na(rotated), 
                        'stationary', 'rotation') 
  
  # Set time point to translation based the information saved by unity
  tra_rot_sta[moving2] <- 'translation'
  
  # Return
  return(tra_rot_sta)
}

# Calculate the state for each time points for each subject and each trial
locomotion_data <- ddply(locomotion_data, c("ppid", "trial"), 
                  mutate, locomotion = what_state(rot_y, moving))

# Calculate the total percentage for each subject
locomotion_per <- ddply(locomotion_data, c("ppid"), summarise,
                        translation = mean(locomotion == "translation"),
                        rotation = mean(locomotion == "rotation"),
                        stationary = mean(locomotion == "stationary"))

# Convert from wide to long format
locomotion_per_long  <- reshape2::melt(locomotion_per, id.vars = c("ppid"))

# Convert proportions to percentage
locomotion_per_long$value <- locomotion_per_long$value * 100

# Write .csv
write.csv(x = locomotion_per_long, file = "other_stuff/source_data_files/locomotion_per_long.csv", 
          quote = FALSE, row.names = FALSE)

# Plot showing percentage of locomotion
locomotion_states_plot <- ggplot(locomotion_per_long, aes(x = variable, y = value)) +
  geom_line(aes(group = ppid), colour = boxplot_line_col) +
  geom_point(colour = boxplot_point_col) +
  geom_boxplot(colour = boxplot_border_col, outlier.shape = NA, width = 0.4,
               fill = base_col[2]) +
  stat_summary(geom = "point", fun = "mean", col = 'white', shape = 24,
               fill = base_col[2], position = position_dodge(width =  0.75), size = 0.8) +
  theme_journal() + 
  #scale_x_discrete(guide = guide_axis(n.dodge = 2)) +
  scale_y_continuous(breaks = c(0, 40, 80)) + 
  coord_cartesian(xlim = c(0.5, 3.5), 
                  ylim = c(0, 80), expand = FALSE) +
  theme(legend.position = "none", plot.margin = unit(c(1, 1, 1, 1), "mm")) +
  labs(title = "", x = "Locomotion", y = "Percentage")

Novelty score over time model

# Load the models
load("fitted_brms_models/SpaNov_m_noveltyScore.Rdata")
# Background: https://www.andrewheiss.com/blog/2022/11/29/conditional-marginal-marginaleffects/
# Set range
x1_range   <- range(m_noveltyScore_run1$data$s_onset_rel)
x1_points  <- seq(from = x1_range[1], to = x1_range[2], length.out = 3)
x2_range   <- range(m_noveltyScore_run2$data$s_onset_rel) 
x2_points  <- seq(from = x2_range[1], to = x2_range[2], length.out = 3)

# Get subject-level predictions
pred1_df <- predictions(m_noveltyScore_run1, newdata = datagrid(s_onset_rel = x1_points, subject = unique), 
                           by = c("s_onset_rel", "subject"))
pred2_df <- predictions(m_noveltyScore_run2, newdata = datagrid(s_onset_rel = x2_points, subject = unique), 
                           by = c("s_onset_rel", "subject"))
pred1_df <- as.data.frame(pred1_df)
pred2_df <- as.data.frame(pred2_df)

# Calculate the minimum for each subject, which is used for colouring
pred1_df <- ddply(pred1_df, c("subject"), mutate, min = min(estimate))
pred2_df <- ddply(pred2_df, c("subject"), mutate, min = min(estimate))

# Create plots
## Run 1
p1 <- ggplot(data = pred1_df, aes(x = s_onset_rel, y = estimate, group = subject, colour = min)) +
  geom_line() +
  scale_color_viridis_c() + theme_journal() +
  theme(legend.position = "none",
        plot.margin = unit(c(1, 2, 1, 2), "mm")) +
  labs(title = "Run 1", x = "Time", y = "Novelty score") +
  scale_x_continuous(breaks = x1_range, labels = c("Start", "End")) +
  scale_y_continuous(breaks = c(-2.0, -0.7,  0.6)) +
  coord_cartesian(xlim = x1_range, ylim = c(-2.03, 0.6),
                  expand = FALSE) 

## Run 2
p2 <- ggplot(data = pred2_df, aes(x = s_onset_rel, y = estimate, group = subject, colour = min)) +
  geom_line() +
  scale_color_viridis_c() + theme_journal() +
  theme(legend.position = "none",
        plot.margin = unit(c(1, 2, 1, 2), "mm")) +
  labs(title = "Run 2", x = "Time", y = "Novelty score") +
  scale_x_continuous(breaks = x2_range, labels = c("Start", "End")) +
  scale_y_continuous(breaks = c(-1.5, -0.5, 0.5)) +
  coord_cartesian(xlim = x2_range, ylim = c(-1.53, 0.5),
                  expand = FALSE) 

# Combine and save
novelty_score_model_plot <- plot_grid(p1, p2, ncol = 2)

Combine to one figure

# Combine
part1 <- plot_grid(total_path_lengths_plot, event_duration_plot, 
                   sector_rayleigh_plot, locomotion_states_plot, align = "h", nrow = 1)
part2 <- plot_grid(NULL, novelty_score_model_plot)
Figure1_combined <- plot_grid(part1, NULL, part2, align = "v", nrow = 3, rel_heights = c(1, 0.2, 1))

# Save as .png
ggsave(Figure1_combined,
       filename = paste0(figurePath, "Fig_1_lower_part.png"), 
       dpi = 1000,
       width = 170,
       height = 90,
       units = "mm")

# Save as .svg
ggsave(Figure1_combined,
       filename = paste0(figurePath, "Fig_1_lower_part.svg"), 
       dpi = 1000,
       width = 170,
       height = 90,
       units = "mm")

Descriptive statistics of exploration

Mean and SD of Raleigh vector length

# Report mean
str1 <- mean_SD_str2(sector_rayleigh$resultant_length, report_type, digits1, rounding_type)
  • Raleigh vector length: M = 0.37 (SD = 0.23)

Calculating the number of visits per participant

# Calculate the number of visited sectors
num_visitedSectors <- ddply(no_control, c("ppid"), summarise, 
                            number = length_uniq(sector))

# Make report string
val  <- num_visitedSectors$number
str1 <- mean_SD_str2(val, report_type, digits1, rounding_type)

Number of visited sector per participant: M = 50 (SD = 2)

Describe two measures on which the novelty score is based

# Load image generated when creating the event files
load("event_tables/images/SpaNov_event_file_gradients.RData")

# Convert list to data frame
data <- as.data.frame(rbindlist(subj_list, idcol = "subject"))

# Add subjects' R number
data$ppid <- subjIDs_R[data$subject]

# Function to match runStartTime
find_runStartTime <- function(ppid, run){
  # Get corresponding to find the run in the trial data
  anonKey <- lookupTable$anonKey[lookupTable$Rnum == ppid[1]]
  
  # Use the anonKey & run to get runStartTime
  bool_index  <- OLM_7T_trial_results$subject == anonKey & 
                 OLM_7T_trial_results$block_num == run
  runStartTime <- OLM_7T_trial_results$runStartTime[bool_index]
  
  return(runStartTime[1])
}

# Use the function to find the correct run start time
data <- ddply(data, c("subject", "ppid", "run"), mutate, 
              runStartTime = find_runStartTime(ppid, run))

# Make time relative to the start of the run. The real onset times will be 
# slightly different but this will not matter for this. 
data$onset_rel <- data$onset - data$runStartTime

# Add run type
data$runType <- "encoding"
data$runType[data$run == 2 | data$run == 4] <- "retrieval"

# Subset to only encoding
data_sub <- data[data$runType == "encoding", ]

# Change run number to match the description in paper. Run 3 is Encoding Run 2
data_sub$run[data_sub$run == 3] <- 2

# Convert run to factor
data_sub$f_run <- as.factor(data_sub$run)


# Calculate descriptive stats for the measures
data_sub_agg1 <- ddply(data_sub, c("subject"), summarise,
                      max_visits = max(visits),
                      min_lastVisit = min(lastVisit, na.rm = TRUE),
                      max_lastVisit = max(lastVisit, na.rm = TRUE),
                      median_lastVisit = median(lastVisit, na.rm = TRUE),
                      mean_lastVisit = mean(lastVisit, na.rm = TRUE))

# Calculate the average duration of events from lastVisits for that exclude NA values
data_sub_agg2 <- ddply(na.omit(data_sub), c("subject"), summarise,
                       mean_duration = mean(duration),
                       median_duration = median(duration),
                       min_duration = min(duration),
                       max_duration = max(duration))

# Create strings for report
val  <- data_sub_agg1$max_visits
str1 <- mean_SD_str2(val, report_type, digits1, rounding_type)
val  <- data_sub_agg1$min_lastVisit
str2 <- mean_SD_str2(val, report_type, digits1, rounding_type, "s")
val  <- data_sub_agg1$max_lastVisit
str3 <- mean_SD_str2(val, report_type, digits1, rounding_type, "s")
val  <- data_sub_agg2$mean_duration
str4 <- mean_SD_str2(val, report_type, digits1, rounding_type, "s")
val  <- data_sub_agg2$min_duration
str5 <- mean_SD_str2(val, report_type, digits1, rounding_type, "s")
val  <- data_sub_agg2$max_duration
str6 <- mean_SD_str2(val, report_type, digits1, rounding_type, "s")
  • Maximum visit per sector for the participants was M = 18 (SD = 1.6) by the time of the second run.
  • The minimum time between visits varied between participants M = 6.7 s (SD = 2.3 s).
  • The maximum time between visits varied between participants M = 1300 s (SD = 310 s).
  • The average duration of the events for the time since varied between participants M = 2.5 s (SD = 0.4 s) with an average minimum of M = 0.6 s (SD = 0.0023 s) and an average maximum of M = 12 s (SD = 5.3 s).

Locomotion during exploration

Report the average values for the three locomotion states

val  <- locomotion_per$translation * 100
str1 <- mean_SD_str2(val, report_type, digits1, rounding_type, "%")
val  <- locomotion_per$rotation  * 100
str2 <- mean_SD_str2(val, report_type, digits1, rounding_type, "%")
val  <- locomotion_per$stationary * 100
str3 <- mean_SD_str2(val, report_type, digits1, rounding_type, "%")
  • Translation: M = 55 % (SD = 8.2 %)
  • Rotation: M = 24 % (SD = 4 %)
  • Stationary: M = 21 % (SD = 6.3 %)

Travelled path length

# Create the report string
str1 <- mean_SD_str2(total_path_lengths$distance, 
                     report_type, digits1, rounding_type, "vm")

The average total path length was M = 3200 vm (SD = 140 vm)

Event duration

# Make report string
val  <- event_df$duration
str1 <- mean_SD_str2(val, report_type, digits1, rounding_type)

Event duration: M = 2.5 (SD = 2.1)

Bayesian hierarchical modelling of novelty score

# Get marginal results
mae_respone_noveltyScore_run1 <- avg_comparisons(m_noveltyScore_run1, variables = "s_onset_rel", 
                                                 type = "response", re_formula = NULL)
mae_respone_noveltyScore_run2 <- avg_comparisons(m_noveltyScore_run2, variables = "s_onset_rel", 
                                                 type = "response", re_formula = NULL)

str1 <- create_str_from_avg_comparisons(mae_respone_noveltyScore_run1, "SDs")
str2 <- create_str_from_avg_comparisons(mae_respone_noveltyScore_run2, "SDs")
  • Average marginal effect in Run 1: -0.48 SDs (95 % CI [-0.55, -0.41])
  • Average marginal effect in Run 2: -0.71 SDs (95 % CI [-0.75, -0.66])

\[ \begin{align*} y & \sim \operatorname{Student}(\nu, \mu, \sigma) \\ \mu & = time + (time | subject) \end{align*} \]

Alternative way to describe: https://rpsychologist.com/r-guide-longitudinal-lme-lmer

Mapping spatial novelty across the hippocampal long axis

Load the results from the linear contrast (Level 1 > Level 2 > Level 3 > …). All xiis from this contrast start with GLM1.

  • Positive values: Higher activity for familiar sectors (c1)
  • Negative values: Higher activity for novel sectors (c2)
# Folder where the images are
ciftiFolder <- "/SpaNov/OLMe_7T_SpaNov_gradient_6lvl_cue-delay_smo4_MSMAll/cope7.feat/stats/vwc/"

# Other parameters
minClusterSize   <- 20

# Choose the files
zMap_file        <- paste0(path2imaging_results2, ciftiFolder, "results_lvl2cope1_dat_ztstat_c1.dscalar.nii")
pMap1_file       <- paste0(path2imaging_results2, ciftiFolder, "results_lvl2cope1_dat_ztstat_cfdrp_c1.dscalar.nii")
pMap2_file       <- paste0(path2imaging_results2, ciftiFolder, "results_lvl2cope1_dat_ztstat_cfdrp_c2.dscalar.nii")
clusterMap1_file <- paste0(path2imaging_results2, ciftiFolder, "results_lvl2cope1_dat_ztstat_cfdrp_c1_clusters.dscalar.nii")
clusterMap2_file <- paste0(path2imaging_results2, ciftiFolder, "results_lvl2cope1_dat_ztstat_cfdrp_c2_clusters.dscalar.nii")
betaMap_file     <- paste0(path2imaging_results2, ciftiFolder, "Y1.dtseries.nii")

# Load all maps as xiis
GLM1_zMap_xii        <- read_cifti(zMap_file, brainstructures = "all")
GLM1_pMap1_xii       <- read_cifti(pMap1_file, brainstructures = "all")
GLM1_pMap2_xii       <- read_cifti(pMap2_file, brainstructures = "all")
GLM1_clusterMap1_xii <- read_cifti(clusterMap1_file, brainstructures = "all")
GLM1_clusterMap2_xii <- read_cifti(clusterMap2_file, brainstructures = "all")
GLM1_betaMap_xii     <- read_cifti(betaMap_file, brainstructures = "all")

# Create masks based on significant effects
HC_familiarity_mask   <- GLM1_pMap1_xii$data$subcort > cutOff & str_detect(GLM1_pMap1_xii$meta$subcort$labels, pattern = "Hippocampus")
HC_novelty_mask       <- GLM1_pMap2_xii$data$subcort > cutOff & str_detect(GLM1_pMap1_xii$meta$subcort$labels, pattern = "Hippocampus")
HC_familiarity_mask_L <- GLM1_pMap1_xii$data$subcort > cutOff & str_detect(GLM1_pMap1_xii$meta$subcort$labels, pattern = "Hippocampus-L")
HC_novelty_mask_L     <- GLM1_pMap2_xii$data$subcort > cutOff & str_detect(GLM1_pMap1_xii$meta$subcort$labels, pattern = "Hippocampus-L")
HC_familiarity_mask_R <- GLM1_pMap1_xii$data$subcort > cutOff & str_detect(GLM1_pMap1_xii$meta$subcort$labels, pattern = "Hippocampus-R")
HC_novelty_mask_R     <- GLM1_pMap2_xii$data$subcort > cutOff & str_detect(GLM1_pMap1_xii$meta$subcort$labels, pattern = "Hippocampus-R")
WB_familiarity_mask   <- GLM1_pMap1_xii > cutOff 
WB_novelty_mask       <- GLM1_pMap2_xii > cutOff 

# Extract & average values
HC_familiarity_values <- GLM1_betaMap_xii$data$subcort[HC_familiarity_mask, ]
HC_familiarity_values <- colMeans(HC_familiarity_values)
HC_novelty_values     <- GLM1_betaMap_xii$data$subcort[HC_novelty_mask, ]
HC_novelty_values     <- colMeans(HC_novelty_values)

HC_familiarity_values_L <- GLM1_betaMap_xii$data$subcort[HC_familiarity_mask_L, ]
HC_familiarity_values_L <- colMeans(HC_familiarity_values_L)

HC_familiarity_values_R <- GLM1_betaMap_xii$data$subcort[HC_familiarity_mask_R, ]
HC_familiarity_values_R <- colMeans(HC_familiarity_values_R)
HC_novelty_values_R     <- GLM1_betaMap_xii$data$subcort[HC_novelty_mask_R, ]
HC_novelty_values_R     <- colMeans(HC_novelty_values_R)

WB_novelty_values <- as.matrix(GLM1_betaMap_xii)
WB_novelty_values <- WB_novelty_values[as.matrix(WB_novelty_mask) == 1, ]
WB_novelty_values <- colMeans(WB_novelty_values)

WB_familiarity_values <- as.matrix(GLM1_betaMap_xii)
WB_familiarity_values <- WB_familiarity_values[as.matrix(WB_familiarity_mask) == 1, ]
WB_familiarity_values <- colMeans(WB_familiarity_values)

# Create cluster tables based on the maps
GLM1_cluster1 <- cifti_cluster_report(zMap_file, 
                            clusterMap1_file, 
                            surfLeft, 
                            surfRight,
                            parcellationFile, 
                            minClusterSize,
                            FALSE)
GLM1_cluster2 <- cifti_cluster_report(zMap_file, 
                            clusterMap2_file, 
                            surfLeft, 
                            surfRight,
                            parcellationFile, 
                            minClusterSize,
                            FALSE)
# Round values to in order to report them
GLM1_cluster1$cluster_values$zValue_max   <- signif(GLM1_cluster1$cluster_values$zValue_max, digits1)
GLM1_cluster1$cluster_values$zValue_min   <- signif(GLM1_cluster1$cluster_values$zValue_min, digits1)
GLM1_cluster1$cluster_values$cluster_mass <- signif(GLM1_cluster1$cluster_values$cluster_mass, digits1)

GLM1_cluster2$cluster_values$zValue_max   <- signif(GLM1_cluster2$cluster_values$zValue_max, digits1)
GLM1_cluster2$cluster_values$zValue_min   <- signif(GLM1_cluster2$cluster_values$zValue_min, digits1)
GLM1_cluster2$cluster_values$cluster_mass <- signif(GLM1_cluster2$cluster_values$cluster_mass, digits1)

# Write tsv files so it's easier to edit
## Positive clusters
### Combine to one table
region_labels  <- GLM1_cluster1$cluster_labels[, 4]
combined_table <- cbind(GLM1_cluster1$cluster_values, region_labels)

### Add places between the region labels
combined_table$region_labels <- str_replace_all(combined_table$region_labels, 
                                                pattern = ",",
                                                replacement = ", ")

### Write as .txt file
write.table(combined_table, file = "figures/SpaNov/tables/GLM1_c1.txt", 
            quote = FALSE,
            row.names = FALSE,
            sep = '\t')

## Negative clusters
### Combine to one table
region_labels  <- GLM1_cluster2$cluster_labels[, 4]
combined_table <- cbind(GLM1_cluster2$cluster_values, region_labels)

### Add places between the region labels
combined_table$region_labels <- str_replace_all(combined_table$region_labels, 
                                                pattern = ",",
                                                replacement = ", ")

### Write as .txt file
write.table(combined_table, file = "figures/SpaNov/tables/GLM1_c2.txt", 
            quote = FALSE,
            row.names = FALSE,
            sep = '\t')

Create unthresholded hippocampal z-map

For unthresholded visualisation of the hippocampal GLM results, we create a version of the z-map that set everything other than the hippocampus to zero.

# (This chunk only needs to run once, so eval is set to FALSE)
# Extreme comparison: Level 1 vs. Level 6
ciftiFolder <- "/SpaNov/OLMe_7T_SpaNov_gradient_6lvl_cue-delay_smo4_MSMAll/cope14.feat/stats/vwc/"
zMap_file1  <- paste0(path2imaging_results2, ciftiFolder, "results_lvl2cope1_dat_ztstat_c1.dscalar.nii")
zMap_xii1   <- read_cifti(zMap_file1, brainstructures = "all", 
                  surfL_fname = surfLeft, surfR_fname = surfRight)

# Linear contrast from Level 1 to Level 6
ciftiFolder <- "/SpaNov/OLMe_7T_SpaNov_gradient_6lvl_cue-delay_smo4_MSMAll/cope7.feat/stats/vwc/"
zMap_file2  <- paste0(path2imaging_results2, ciftiFolder, "results_lvl2cope1_dat_ztstat_c1.dscalar.nii")
zMap_xii2   <- read_cifti(zMap_file2, brainstructures = "all", 
                  surfL_fname = surfLeft, surfR_fname = surfRight)

# Create mask for everything other than right and left hippocampus
currentMask <- zMap_xii1$meta$subcort$labels != "Hippocampus-R" & 
               zMap_xii1$meta$subcort$labels != "Hippocampus-L"

# Create new variables
zMap_xii1_HC <- zMap_xii1
zMap_xii2_HC <- zMap_xii2

# Set everything within this mask to zero
zMap_xii1_HC$data$subcort[currentMask, ] <- 0
zMap_xii2_HC$data$subcort[currentMask, ] <- 0

# Also set cortical values to zero
zMap_xii1_HC$data$cortex_left  <- matrix(as.integer(0), nrow = 29696, ncol = 1)
zMap_xii1_HC$data$cortex_right <- matrix(as.integer(0), nrow = 29716, ncol = 1)
zMap_xii2_HC$data$cortex_left  <- matrix(as.integer(0), nrow = 29696, ncol = 1)
zMap_xii2_HC$data$cortex_right <- matrix(as.integer(0), nrow = 29716, ncol = 1)

# Create new names
new_zMap_file1 <- str_replace(zMap_file1, pattern = ".dscalar.nii", 
                              replacement = "_onlyHC.dscalar.nii")
new_zMap_file2 <- str_replace(zMap_file2, pattern = ".dscalar.nii", 
                              replacement = "_onlyHC.dscalar.nii")

# Write new cifti files
write_cifti(xifti = zMap_xii1_HC, cifti_fname = new_zMap_file1)
write_cifti(xifti = zMap_xii2_HC, cifti_fname = new_zMap_file2)

Spatial novelty-gradient in hippocampus

Calculate GS & tSNR in hippocampus

# Load tSNR and GS values
load("data/ignore_fMRI_version1/tSNR_GS_maps.RData")

# Add file ID
GS_tSNR_df$rowID <- 1:nrow(GS_tSNR_df)

# Subset to subjects included in the analysis
GS_tSNR_df <- GS_tSNR_df[GS_tSNR_df$GS_subject %in% subjIDs_R, ]

# Calculate average GS & tSNR maps
GS_merged   <- merge_xifti(xifti_list = GS_list[GS_tSNR_df$rowID])
GS_avg      <- apply_xifti(GS_merged, margin = 1, mean)
tSNR_merged <- merge_xifti(xifti_list = tSNR_list[GS_tSNR_df$rowID])
tSNR_avg    <- apply_xifti(tSNR_merged, margin = 1, mean)

# Get MNI coordinates
HC_data_L  <- MNI_coord[str_starts(MNI_coord$region, pattern = "Hippocampus-L"), ]
HC_data_R  <- MNI_coord[str_starts(MNI_coord$region, pattern = "Hippocampus-R"), ]
HC_data    <- rbind(HC_data_L, HC_data_R)

# Get bool index 
HC_index_L <- GS_list[[1]]$meta$subcort$labels == "Hippocampus-L"
HC_index_R <- GS_list[[1]]$meta$subcort$labels == "Hippocampus-R"

# Results list
GS_HC_list   <- list()
tSNR_HC_list <- list()

# Extract hippocampus values by loop over all subjects
for(i in 1:length(subjIDs_R)){
  # GS
  ## Extract values from Run 1 & 2
  run1_xii <- GS_list[[GS_tSNR_df$rowID[GS_tSNR_df$GS_run == "RUN1" & GS_tSNR_df$GS_subject == subjIDs_R[i]]]]
  run2_xii <- GS_list[[GS_tSNR_df$rowID[GS_tSNR_df$GS_run == "RUN2" & GS_tSNR_df$GS_subject == subjIDs_R[i]]]]
  run1_HC  <- c(run1_xii$data$subcort[HC_index_L], run1_xii$data$subcort[HC_index_R])
  run2_HC  <- c(run2_xii$data$subcort[HC_index_L], run2_xii$data$subcort[HC_index_R])
  
  ## Average across runs
  avg_HC   <- (run1_HC + run2_HC)/2
  
  ## Add to data.frame
  temp_df         <- HC_data
  temp_df$subject <- subjIDs_R[i]
  temp_df$GS      <- avg_HC
  GS_HC_list[[i]] <- temp_df
  
  # GS
  ## Extract values from Run 1 & 2
  run1_xii <- tSNR_list[[GS_tSNR_df$rowID[GS_tSNR_df$GS_run == "RUN1" & GS_tSNR_df$GS_subject == subjIDs_R[i]]]]
  run2_xii <- tSNR_list[[GS_tSNR_df$rowID[GS_tSNR_df$GS_run == "RUN2" & GS_tSNR_df$GS_subject == subjIDs_R[i]]]]
  run1_HC  <- c(run1_xii$data$subcort[HC_index_L], run1_xii$data$subcort[HC_index_R])
  run2_HC  <- c(run2_xii$data$subcort[HC_index_L], run2_xii$data$subcort[HC_index_R])
  
  ## Average across runs
  avg_HC   <- (run1_HC + run2_HC)/2
  
  ## Add to data.frame
  temp_df           <- HC_data
  temp_df$subject   <- subjIDs_R[i]
  temp_df$tSNR      <- avg_HC
  tSNR_HC_list[[i]] <- temp_df
}

# Convert to data frame
GS_HC_df   <- rbindlist(GS_HC_list)
tSNR_HC_df <- rbindlist(tSNR_HC_list)
GS_HC_df   <- as.data.frame(GS_HC_df)
tSNR_HC_df <- as.data.frame(tSNR_HC_df)

# Simple A vs. P comparison
## Divide into A & P
GS_HC_df$AP   <- ifelse(GS_HC_df$y > -21, "Anterior", "Posterior")
tSNR_HC_df$AP <- ifelse(tSNR_HC_df$y > -21, "Anterior", "Posterior")

## Calculate average for each participant
GS_HC_df_AP   <- ddply(GS_HC_df, c("subject", "AP"), summarise, GS = mean(GS))
tSNR_HC_df_AP <- ddply(tSNR_HC_df, c("subject", "AP"), summarise, tSNR = mean(tSNR))

## Run tests
### GS
GS_A <- GS_HC_df_AP$GS[GS_HC_df_AP$AP == "Anterior"]
GS_P <- GS_HC_df_AP$GS[GS_HC_df_AP$AP == "Posterior"]
diff <- GS_A - GS_P
GS_d    <- round(mean(diff)/sd(diff), 2)
GS_p    <- t.test(diff)$p.value
GS_BF10 <- signif(reportBF(ttestBF(diff)), 3)

### tSNR
tSNR_A <- tSNR_HC_df_AP$tSNR[tSNR_HC_df_AP$AP == "Anterior"]
tSNR_P <- tSNR_HC_df_AP$tSNR[tSNR_HC_df_AP$AP == "Posterior"]
diff      <- tSNR_A - tSNR_P
tSNR_d    <- round(mean(diff)/sd(diff), 2)
tSNR_p    <- t.test(diff)$p.value
tSNR_BF10 <- signif(reportBF(ttestBF(diff)), 3)

# Calculate percentile of averge hippocampus tSNR
tSNR_A_percentile <- round(mean(as.matrix(tSNR_avg) < mean(tSNR_A)), 2)
tSNR_P_percentile <- round(mean(as.matrix(tSNR_avg) < mean(tSNR_P)), 2)

# Get values from the XIFTIs
HC_tSNR_avg <- c(tSNR_avg$data$subcort[HC_index_L], tSNR_avg$data$subcort[HC_index_R])
HC_GS_avg   <- c(GS_avg$data$subcort[HC_index_L], GS_avg$data$subcort[HC_index_R])

## Create a data frame
HC_tSNR_GS_df   <- data.frame(tSNR = HC_tSNR_avg, GS = HC_GS_avg, 
                              region = rep(c("Hippocampus-L", "Hippocampus-R"), 
                                          times = c(sum(HC_index_L), sum(HC_index_R))))

Permutation analyses

# Load hippocampal gradient data
load("intermediate_data/SpaNov_gradient_data_cue-delay.RData")

# Add tSNR & GS
HC_data$tSNR <- HC_tSNR_GS_df$tSNR
HC_data$GS   <- HC_tSNR_GS_df$GS
# In 5b1e712, carefully check that the order HC_data and HC_tSNR_GS_df is 
# the same. Specifically I calculated the position values for HC_tSNR_GS_df
# and then checked if mean(HC_data$position == HC_tSNR_GS_df$position) == 1

# Average across the AP position
HC_data_agg_pos <- ddply(HC_data, c("position", "Hemisphere"), summarise, 
                         n = length(min), min = mean(min), max = mean(max), 
                         tSNR = mean(tSNR), GS = mean(GS))

# Average across position and hemisphere
HC_data_agg_pos2 <- ddply(HC_data, c("position"), summarise, 
                          n = length(min), min = mean(min), max = mean(max),
                          tSNR = mean(tSNR), GS = mean(GS))

# Labels
conN       <- 6
novLabels  <- paste0('Level ', 1:conN)
permutation_analysis <- function(data, lm_formula, nIter, colName, imageName){
  # Initialise results list
  results <- list()
  
  # Select correct column for analysis and create new data frame
  data$val     <- data[, colName]
  data2shuffle <- data
  
  # Calculate empirical values and save to list
  results$lm <- lm(lm_formula, data = data)
  numCoef    <- length(results$lm$coefficients) - 1 # Ignoring the intercept
  
  # Start cluster
  my.cluster <- parallel::makeCluster(detectCores() - 2, type = "PSOCK")
  
  #register it to be used by %dopar%
  doParallel::registerDoParallel(cl = my.cluster)
  
  # Run parallel loop
  permuted_values <- foreach(i = 1:nIter, .combine = 'c', .packages = 'plyr') %dopar% {
    data2shuffle$val <- sample(data$val)
    
    # Fit model
    temp_lm <- lm(lm_formula, data = data2shuffle)
    
    # Add values 
    temp_est <- as.data.frame(matrix(as.numeric(temp_lm$coefficients)[-1], ncol = numCoef))
    names(temp_est) <- names(results$lm$coefficients)[-1]
    list(temp_est)
  }
  
  # Stop cluster again
  parallel::stopCluster(cl = my.cluster) 
  
  # Add to results
  results$permuted_values <- as.data.frame(rbindlist(permuted_values))
  
  # Save to disk
  if(!missing(imageName)){
    save(results, file = paste0("intermediate_data/", imageName))
  }
  
  # Return value
  return(results)
}
Minimum
# Rename to be unique
grad_HC1 <- HC_data_agg_pos

# Check if the code needs to be run again. This is done by checking the md5 has 
# for the data frame that is used in the calculation if it is different from the 
# hash sum saved in md5_hash_table.csv, it is re-run. This is to avoid having to 
# re-run everything each time I work on the data. 
runCodeAgain1 <- check_if_md5_hash_changed(grad_HC1, hash_table_name = "SpaNov_md5_hash_table.csv")
Model with interaction with hemisphere
# Seed
set.seed(19911225)

# Other input
lm_formula    <- "val ~ position * Hemisphere + tSNR + GS"
nIter         <- 100000
colName       <- "min"
imageName     <- "SpaNov_permut_HC_grad_analysis1_cue-delay.RData"

# Run if necessary
if(runCodeAgain1){
  grad_min_permut1 <- permutation_analysis(grad_HC1, lm_formula, 
                                       nIter, colName, imageName)
} else {
  load(paste0("intermediate_data/", imageName))
  grad_min_permut1 <- results
}

# Select values for plotting
dist                <- grad_min_permut1$permuted_values[,5]
critVal             <- grad_min_permut1$lm$coefficients[6]
grad_min_permut1_p5 <- pValue_from_nullDist(critVal, dist, "two.sided")
Left hippocampus
# Seed
set.seed(19911225)

# Other input
lm_formula    <- "val ~ position + tSNR + GS"
nIter         <- 100000
colName       <- "min"
imageName     <- "SpaNov_permut_HC_grad_analysis1_cue-delay_L.RData"

# Run if necessary
if(runCodeAgain1){
  grad_min_permut1_L <- permutation_analysis(grad_HC1[grad_HC1$Hemisphere == 'left', ], lm_formula, 
                                       nIter, colName, imageName)
} else {
  load(paste0("intermediate_data/", imageName))
  grad_min_permut1_L <- results
}

# Select values for plotting
dist                  <- grad_min_permut1_L$permuted_values[,1]
critVal               <- grad_min_permut1_L$lm$coefficients[2]
grad_min_permut1_L_p1 <- pValue_from_nullDist(critVal, dist, "two.sided")
Right hippocampus
# Seed
set.seed(19911225)

# Other input
lm_formula    <- "val ~ position + tSNR + GS"
nIter         <- 100000
colName       <- "min"
imageName     <- "SpaNov_permut_HC_grad_analysis1_cue-delay_R.RData"

# Run if necessary
if(runCodeAgain1){
  grad_min_permut1_R <- permutation_analysis(grad_HC1[grad_HC1$Hemisphere == 'right', ], lm_formula, 
                                       nIter, colName, imageName)
} else {
  load(paste0("intermediate_data/", imageName))
  grad_min_permut1_R <- results
}

# Select values for plotting
dist                  <- grad_min_permut1_R$permuted_values[,1]
critVal               <- grad_min_permut1_R$lm$coefficients[2]
grad_min_permut1_R_p1 <- pValue_from_nullDist(critVal, dist, "two.sided")
Maximum
Model with interaction with hemisphere
# Seed
set.seed(19911225)

# Other input
lm_formula    <- "val ~ position * Hemisphere + tSNR + GS"
nIter         <- 100000
colName       <- "max"
imageName    <- "SpaNov_permut_HC_grad_analysis3_cue-delay.RData"

# Run if necessary
if(runCodeAgain1){
  grad_max_permut1 <- permutation_analysis(grad_HC1, lm_formula, 
                                       nIter, colName, imageName)
} else {
  load(paste0("intermediate_data/", imageName))
  grad_max_permut1 <- results
}

# Select values for plotting
dist                <- grad_max_permut1$permuted_values[,5]
critVal             <- grad_max_permut1$lm$coefficients[6]
grad_max_permut1_p5 <- pValue_from_nullDist(critVal, dist, "two.sided")
Average across hemispheres
# Rename to be unique
grad_HC2 <- HC_data_agg_pos2

# Check if the code needs to be run again. This is done by checking the md5 has 
# for the data frame that is used in the calculation if it is different from the 
# hash sum saved in md5_hash_table.csv, it is re-run. This is to avoid having to 
# re-run everything each time I work on the data. 
runCodeAgain2 <- check_if_md5_hash_changed(grad_HC2, hash_table_name = "SpaNov_md5_hash_table.csv")

# Seed
set.seed(20131)

# Other input
lm_formula    <- "val ~ position + tSNR + GS"
nIter         <- 100000
colName       <- "max"
imageName    <- "SpaNov_permut_HC_grad_analysis4_cue-delay.RData"

# Run if necessary
if(runCodeAgain2){
  grad_max_permut1_avg <- permutation_analysis(grad_HC2, lm_formula, 
                                       nIter, colName, imageName)
} else {
  load(paste0("intermediate_data/", imageName))
  grad_max_permut1_avg <- results
}

# Select values for plotting
dist              <- grad_max_permut1_avg$permuted_values[,1]
critVal           <- grad_max_permut1_avg$lm$coefficients[2]
grad_max_permut1_avg_p1 <- pValue_from_nullDist(critVal, dist, "two.sided")

Report stat of permutation analysis

# Calculate effect sizes
grad_min_permut1_eff     <- eta_squared(car::Anova(grad_min_permut1$lm, type = 2))
grad_min_permut1_L_eff   <- eta_squared(car::Anova(grad_min_permut1_L$lm, type = 2))
grad_min_permut1_R_eff   <- eta_squared(car::Anova(grad_min_permut1_R$lm, type = 2))
grad_max_permut1_R_eff   <- eta_squared(car::Anova(grad_max_permut1$lm, type = 2))
grad_max_permut1_avg_eff <- eta_squared(car::Anova(grad_max_permut1_avg$lm, type = 2))

# Create data frame
## Create variables
tmp_formula   <- c("min ~ position * Hemisphere + tSNR + GS", "min (left) ~ position + tSNR + GS",
                   "min (right) ~ position + tSNR + GS", "max ~ position * Hemisphere + tSNR + GS",
                   "max (avg) ~ position + tSNR + GS")
tmp_coef_name <- c("interaction", "position", "position", "interaction", "position")
tmp_coef_val  <- c(grad_min_permut1$lm$coefficients[6], grad_min_permut1_L$lm$coefficients[2],
                   grad_min_permut1_R$lm$coefficients[2], grad_max_permut1$lm$coefficients[6],
                   grad_max_permut1_avg$lm$coefficients[2])
tmp_coef_es   <- c(grad_min_permut1_eff$Eta2_partial[5], grad_min_permut1_L_eff$Eta2_partial[1],
                   grad_min_permut1_R_eff$Eta2_partial[1], grad_max_permut1_R_eff$Eta2_partial[5],
                   grad_max_permut1_avg_eff$Eta2_partial[1])
tmp_coef_p    <- c(grad_min_permut1_p5, grad_min_permut1_L_p1, grad_min_permut1_R_p1, 
                   grad_max_permut1_p5, grad_max_permut1_avg_p1)

## Convert to numeric values and flip
tmp_coef_val <- -as.numeric(tmp_coef_val)

## Add to one data frame
tab_df <- data.frame(Formula = tmp_formula, Coefficient = tmp_coef_name,
                     beta = tmp_coef_val, eta_squared = tmp_coef_es, p = tmp_coef_p)

## Round columns
tab_df$beta        <- round(tab_df$beta, 3)
tab_df$eta_squared <- round(tab_df$eta_squared, 3)

## Create p-values by looping over all values
for(i in seq_along(tab_df$p)){
  tab_df$p[i] <-paste("p", pValue(as.numeric(tab_df$p[i])))
}

# Show table
kable(tab_df)
Formula Coefficient beta eta_squared p
min ~ position * Hemisphere + tSNR + GS interaction 0.022 0.030 p = .021
min (left) ~ position + tSNR + GS position 0.021 0.053 p = .021
min (right) ~ position + tSNR + GS position 0.039 0.154 p < .001
max ~ position * Hemisphere + tSNR + GS interaction 0.013 0.011 p = .194
max (avg) ~ position + tSNR + GS position 0.031 0.104 p < .001

Indidiviudal slopes

# Load individual slopes
load("data/ignore_fMRI_version1/extracted_values/SpavNov_gradient_subject-slopes_20250721_073753.RData")

# Add significance
subjSlopes$sig <- subjSlopes$`Pr(>|t|)` < .05

# Individual slopes from gradient analysis
gradient_subject_slopes_hemisphere     <- -subjSlopes$Estimate[subjSlopes$type == "position" & subjSlopes$effect == "position:Hemisphereright"]
gradient_subject_slopes_position_L     <- -subjSlopes$Estimate[subjSlopes$type == "position_L" & subjSlopes$effect == "position"]
gradient_subject_slopes_position_R     <- -subjSlopes$Estimate[subjSlopes$type == "position_R" & subjSlopes$effect == "position"]
gradient_subject_slopes_position_L_sig <- subjSlopes$sig[subjSlopes$type == "position_L" & subjSlopes$effect == "position"]
gradient_subject_slopes_position_R_sig <- subjSlopes$sig[subjSlopes$type == "position_R" & subjSlopes$effect == "position"]

# Effect of hemispheres
x      <- gradient_subject_slopes_hemisphere
d1     <- signif(mean(x)/sd(x), digits1)
str1   <- mean_SD_str2(x, report_type, digits1, rounding_type)
BF10_1 <- reportBF(ttestBF(x))

# Effect in left hemisphere
x      <- gradient_subject_slopes_position_L
d3     <- signif(mean(x)/sd(x), digits1)
str3   <- mean_SD_str2(x, report_type, digits1, rounding_type)
BF10_3 <- reportBF(ttestBF(x))

# Effect in right hemisphere
x      <- gradient_subject_slopes_position_R
d4     <- signif(mean(x)/sd(x), digits1)
str4   <- mean_SD_str2(x, report_type, digits1, rounding_type)
BF10_4 <- reportBF(ttestBF(x))

In the right hippocampus, individual slopes were different from zero,M = 0.0093 (SD = 0.017), BF10 = 148.74, d = 0.54, while they was weak evidence for that null hypothesis in the left hippocampus,M = 0.0037 (SD = 0.017), BF10 = 0.49, d = 0.21. However, evidence for a difference between hemisphere in terms of individual slopes was inconclusive,M = 0.0056 (SD = 0.017), BF10 = 2.04, d = 0.32.

Figure 2

# Create new data frame just for plotting
HC_data_plotting <- HC_data_agg_pos
HC_data_plotting$Hemisphere2 <- ifelse(HC_data_plotting$Hemisphere == "right",
                                       "Right hippocampus", "Left hippocampus")

# Write ,csv file
write.csv(x = HC_data_plotting, file = "other_stuff/source_data_files/HC_scatter_plot.csv", 
          quote = FALSE, row.names = FALSE)

# Create Figure 2
## Create minimum scatter plot for both hemispheres
scatter_min <- ggplot(HC_data_plotting, aes(x = -position, y = min)) + 
  facet_grid(~Hemisphere2, scales = "free_x") + 
  geom_point(aes(colour = min)) +
  geom_smooth(method = "lm", formula = y ~ x, colour = "black") +
  labs(x = "Distance to most anterior part (mm)", y = "Novelty\npreference") + 
  scale_x_continuous(breaks =  seq(from = -45, to = 0, by = 15), labels = c("45", "30", "15", "0")) +
  scale_y_continuous(breaks = 1:6, labels = paste("Lvl", 1:6), limits = c(1, 6)) +
  scale_colour_viridis_c(option = "H", limits = c(1, 6), direction = -1) +
  theme_journal() +
  theme(strip.background = element_rect(color="white", fill="white"),
        axis.title.y = element_text(margin = margin(t = 0, r = 8, b = 0, l = 0)),
        panel.spacing.x = unit(7.5, "mm")) +
  theme(legend.position = "none")

## Create plots for individual slopes
### Create data frame for plotting
individual_slope_df <- data.frame(Effect = rep(c("Right hippocampus", "Left hippocampus", "Interaction"), each = 56),
                                  Coefficient = c(gradient_subject_slopes_position_R,
                                                  gradient_subject_slopes_position_L,
                                                  gradient_subject_slopes_hemisphere))

# Write ,csv file
write.csv(x = individual_slope_df, file = "other_stuff/source_data_files/HC_individual_slopes.csv", 
          quote = FALSE, row.names = FALSE)

### Create plot
individual_slopes_plot <- ggplot(individual_slope_df, aes(x = Effect, y = Coefficient)) +
  geom_hline(yintercept = 0, linetype = 2) +
  geom_jitter(colour = boxplot_point_col, width = 0.25) +
  geom_boxplot(colour = boxplot_border_col, alpha = 0.7, outlier.shape = NA, width = 0.4,
               fill = base_col[2]) +
  stat_summary(geom = "point", fun = "mean", col = 'white', shape = 24,
               fill = base_col[2], position = position_dodge(width =  0.75), size = 0.8) +
  theme_journal() + 
  theme(axis.title.y = element_text(margin = margin(t = 0, r = 4.5, b = 0, l = 0))) +
  labs(title = "", x = "", y = "Individual\nregression coefficients")

## Combine both plots
p <- plot_grid(scatter_min, individual_slopes_plot, ncol = 1, align = "hv", axis = "tblr")

## Save as .png
ggsave(p, filename = paste0(figurePath, "Fig_2c_scattter_and_boxplot.png"),
       dpi = dpi,  width = 99, height = 83.5, units = "mm")

## Save as .svg
ggsave(p, filename = paste0(figurePath, "Fig_2c_scattter_and_boxplot.svg"),
       dpi = dpi,  width = 99, height = 83.5, units = "mm")

Extension of spatial novelty gradient to the posterior medial cortex

Figure 3

Yeo 7 parcellation ridge plot

# Get the names of the networks (-1 to remove the first row that just contains ???)
Yeo7_labels        <- Yeo7_xii$meta$cifti$labels$parcels[-1, ]
Yeo7_labels$Region <- row.names(Yeo7_labels)
row.names(Yeo7_labels)   <- NULL

# Convert the RGB values to hexcodes
Yeo7_labels$hexcol <- rgb(Yeo7_labels$Red, Yeo7_labels$Green, Yeo7_labels$Blue, 
                          maxColorValue = 1)

# Create basic data frame
Yeo_zValues_df <- data.frame(zValue = c(GLM1_zMap_xii$data$cortex_left, GLM1_zMap_xii$data$cortex_right),
                             Yeo7   = c(Yeo7_xii$data$cortex_left, Yeo7_xii$data$cortex_right))

# Add network label
Yeo_zValues_df$net_label <- NA
for(i in 1:nrow(Yeo7_labels)){
  # Select all rows in Yeo_zValues_df that have this key
  bool_index <- Yeo_zValues_df$Yeo7 == Yeo7_labels$Key[i]
  Yeo_zValues_df$net_label[bool_index] <- Yeo7_labels$Region[i]
}

# Find lowest z value that's significant
## Get the p-values and the z-values
GLM1_pValues1 <- get_all_points_from_xifti(GLM1_pMap1_xii)
GLM1_pValues2 <- get_all_points_from_xifti(GLM1_pMap2_xii)
GLM1_zValues  <- get_all_points_from_xifti(GLM1_zMap_xii)

## Subset to only significant values
GLM1_zValues1 <- GLM1_zValues[GLM1_pValues1 > cutOff]
GLM1_zValues2 <- GLM1_zValues[GLM1_pValues2 > cutOff]

# Exclude vertices that are not included in Yeo 7
Yeo_zValues_df_sub <- Yeo_zValues_df[Yeo_zValues_df$Yeo7 != 0, ]

# Order according to the mean of the network
Yeo_zValues_df_agg <- ddply(Yeo_zValues_df_sub, c("net_label"), summarise, 
                            avg_z = mean(zValue))

## First order alphabetically so we can apply the same order to GLM1_Yeo7
Yeo_zValues_df_agg <- Yeo_zValues_df_agg[order(as.character(Yeo_zValues_df_agg$net_label)), ]

## Now order based on the mean
mean_order         <- order(Yeo_zValues_df_agg$avg_z)
Yeo_zValues_df_agg <- Yeo_zValues_df_agg[mean_order, ]

# Order the colours in the same way
Yeo7_labels <- Yeo7_labels[order(as.character(Yeo7_labels$Region)), ]
Yeo7_labels <- Yeo7_labels[mean_order, ]

# Add the full names
Yeo7_labels$fullNames <- find_values_thru_matching(Yeo7_abbr, Yeo7_fullNames, Yeo7_labels$Region)

# Make net_label a factor with the correct order
Yeo_zValues_df_sub$net_label <- factor(Yeo_zValues_df_sub$net_label, 
                                       levels = Yeo_zValues_df_agg$net_label,
                                       labels = Yeo7_labels$fullNames,
                                       ordered = TRUE)

# Create plot
p1 <- ggplot(Yeo_zValues_df_sub, aes(x = zValue, y = net_label, fill = net_label)) + 
  geom_density_ridges(linewidth = 0.3) +
  scale_fill_manual(values = Yeo7_labels$hexcol) +
  geom_vline(xintercept = max(GLM1_zValues2), linetype = 2, linewidth = 0.3) +
  geom_vline(xintercept = min(GLM1_zValues1), linetype = 2, linewidth = 0.3) +
  theme_journal() + 
  coord_cartesian(ylim = c(0.5, 8.5)) +
  theme(legend.position = "none", plot.margin = unit(c(1,1,1,1), "mm")) +
  labs(x = "z-statistic", y = "")

# Save as .png
ggsave(p1,
       filename = paste0(figurePath, "Fig_3a_Yeo7_ridge.png"), 
       dpi = dpi,
       width = 60,
       height = 30,
       units = "mm")

# Save as .svg
ggsave(p1,
       filename = paste0(figurePath, "Fig_3a_Yeo7_ridge.svg"), 
       dpi = dpi,
       width = 60,
       height = 30,
       units = "mm")

Margulies’ connectivity gradients

# Load all 10 gradients from brainstat
# Source: https://brainstat.readthedocs.io/en/master/index.html
nGrad          <- 10
prefix         <- "data/ignore_fMRI_version1/brainstat/Gradient_"

# Get the medial walls from donour
mwm_L <- Yeo7_xii$meta$cortex$medial_wall_mask$left
mwm_R <- Yeo7_xii$meta$cortex$medial_wall_mask$right

# Create data frame with 
Margulies_grad        <- as.data.frame(matrix(NA, nrow = 59412, ncol = 10))
names(Margulies_grad) <- paste0("Grad_", 1:10)

# Loop through all gradients
for(i in 1:nGrad){
  # Create a new xifti by reading the gifti files
  tmp_xii <- read_xifti2(cortexL = paste0(prefix, i, "_L.shape.gii"), 
                                       surfL = surfLeft,
                                       mwall_values = c(0),
                                       cortexR = paste0(prefix, i, "_R.shape.gii"),
                                       surfR = surfRight)
  
  # Use the medial wall values from donour and apply them to the loaded file. 
  tmp_xii$data$cortex_left[!mwm_L]  <- NA
  tmp_xii$data$cortex_right[!mwm_R] <- NA
  tmp_xii <- move_to_mwall(tmp_xii, values = c(NA))
  
  # Add to prepared data frame
  Margulies_grad[, i] <- c(tmp_xii$data$cortex_left, tmp_xii$data$cortex_right)
}

# Add Yeo 7 to the data frame
Margulies_grad$Yeo7_name   <- Yeo_zValues_df$net_label

# Remove the missing values
Margulies_grad_sub <- Margulies_grad[!is.na(Margulies_grad$Yeo7_name), ]

# Plot 
p1 <- ggplot(Margulies_grad_sub, aes(x = Grad_2, y = Grad_1, colour = Yeo7_name)) + 
  geom_point(alpha = 0.2, size = 0.3, shape = 16) +
  theme_void() + 
  theme(legend.position = "none") +
  scale_colour_manual(values = Yeo7_colours)

# Save file
ggsave(p1,
       filename = paste0(figurePath, "Fig_3a_Yeo7_Margulies_space.png"), 
       dpi = dpi,
       width = 30,
       height = 30,
       units = "mm")

# Significant vertices for GLM1
Margulies_grad$GLM_1_direction <- NA
Margulies_grad$GLM_1_direction[c(GLM1_pMap1_xii$data$cortex_left, GLM1_pMap1_xii$data$cortex_right) > cutOff]  <- "Familiarity"
Margulies_grad$GLM_1_direction[c(GLM1_pMap2_xii$data$cortex_left, GLM1_pMap2_xii$data$cortex_right) > cutOff]  <- "Novelty"

# Create subset of only significant vertices
Margulies_grad_sub <- Margulies_grad[!is.na(Margulies_grad$GLM_1_direction), ]

# Plot 
p1 <- ggplot(Margulies_grad_sub, aes(x = Grad_2, y = Grad_1, colour = GLM_1_direction)) + 
  geom_point(alpha = 0.2, size = 0.3, shape = 16) +
  theme_void() + 
  scale_colour_manual(values = cool_warm_colours) + 
  theme(legend.position = "none") 

# Save file
ggsave(p1,
       filename = paste0(figurePath, "Fig_3a_NovFam_Margulies_space.png"), 
       dpi = dpi,
       width = 70,
       height = 70,
       units = "mm")
# Splitting the dataset into the Training set and Test set 
# https://www.geeksforgeeks.org/classifying-data-using-support-vector-machinessvms-in-r/
svm_df <-  data.frame(Grad_1 = Margulies_grad_sub$Grad_1,
                      Grad_2 = Margulies_grad_sub$Grad_2,
                      Modulation = ifelse(Margulies_grad_sub$GLM_1_direction == "Novelty", 1, 0))
# Seed for this analysis
set.seed(123) 

# Split the data into training and test
split        <- sample.split(svm_df$Modulation, SplitRatio = 0.75) 
training_set <- subset(svm_df, split == TRUE) 
test_set     <- subset(svm_df, split == FALSE) 

# Feature scaling 
training_set[-3] <- scale(training_set[-3]) 
test_set[-3]     <- scale(test_set[-3]) 

# Fitting SVM to the Training set 
classifier <- svm(formula = Modulation ~ ., 
                 data = training_set, 
                 type = 'C-classification', 
                 kernel = 'radial') 

# Remove NA
test_set <- na.omit(test_set)

# Predicting the Test set results
y_pred   <- predict(classifier, newdata = test_set) 
cm       <- table(test_set[, 3], y_pred)/ sum(table(test_set[, 3], y_pred))
cm       <- round(cm *100, 2)
accuracy <- mean(test_set[, 3] == y_pred)

The SVM-classification accuracy is 95.17% for significant novelty vs. familiarity vertices in Margulies’ et al. (2015) gradient space.

Resting-state connectivity between spatial novelty gradients in hippocampus and the posterior parietal cortex

Figure 4

############ connectivity_schematic
# Create hippocampus gradient bar
SC_gradient_xii    <- read_cifti("cifti_results/SpaNov_gradient_wholebrain_min_cue-delay.dlabel.nii", brainstructures = c("subcort"))

# Create hippocampus only gradient map and write as CIFTI
HC_gradient_xii <- SC_gradient_xii
HC_gradient_xii$data$subcort[!str_detect(SC_gradient_xii$meta$subcort$labels, pattern = "Hippocampus"), 1] <- as.integer(0)
#write_cifti(HC_gradient_xii, cifti_fname = "cifti_results/SpaNov_gradient_HC_min_cue-delay.dlabel.nii")

############# Model results
# Add Marginal results
load("fitted_brms_models/SpaNov_m_conn2.Rdata")

# Rename model & set variable of interest
m   <- m_conn2
voi <- "diagonal"

# Plot means of the conditions
## Use prediction
pred_means <- predictions(m, newdata = datagrid(diagonal = unique, f_gradient_level_cortex = unique, 
                                                f_gradient_level_HC = unique, subject = unique), 
                           by = voi, re_formula = NULL, type = "response")
pred_means <- as.data.frame(pred_means) 

## Get posterior draws
pred_means_draws <-  get_draws(pred_means)

# Plot means for each participant
pred_means_subject <- avg_comparisons(m, variables = voi, by = c("subject"))
pred_means_subject <- as.data.frame(pred_means_subject) 

## Get posterior draws
pred_means_subject_draws <- get_draws(pred_means_subject)
pred_means_subject_draws <- as.data.frame(pred_means_subject_draws)

# Order by mean difference
subject_ordered <- pred_means_subject$subject[order(pred_means_subject$estimate, decreasing = TRUE)]
pred_means_subject_draws$f_subject <- factor(pred_means_subject_draws$subject,
                                                 levels = subject_ordered, ordered = TRUE)

# Plot
# Width: a vector of probabilities to use that determine the widths of the resulting intervals 
# Here that to 95%
pred_means_subject_draws_p <- ggplot(pred_means_subject_draws, aes(y = f_subject, x = -draw, colour = f_subject)) +
  geom_vline(xintercept = 0, linetype = "dashed") + 
  stat_pointinterval(.width = 0.95, linewidth = 0.25, point_size = 0.05) +
  labs(x = "Within novelty - between\nnovelty connectivity", y = "Individuals") +
  theme_journal() +
  theme(plot.title = element_text(hjust = 0.5),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        plot.subtitle = element_text(hjust = 0.5),
        legend.position = "none")

# Average difference between the conditions
## Calculate marginal average effect is the contrast 0 vs. 1
mae1_respone <- avg_comparisons(m, variables = voi, type = "response", re_formula = NULL)

## Get posterior draws
mae_response_draws <-  get_draws(mae1_respone)

## Flipping sign because it makes more sense
mae_response_draws$draw <- -mae_response_draws$draw

## Create plot
p1_mae_draws <- ggplot(mae_response_draws, aes(x = draw)) +
  geom_vline(xintercept = 0, linetype = "dashed") + 
  stat_halfeye(fill = base_col[2], point_size = 1.5) +
  labs(x = "Within novelty - between\nnovelty connectivity", y = "Density") +
  theme_journal() +
  theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5)) 

# Save as .png
ggsave(plot_grid(p1_mae_draws, pred_means_subject_draws_p, align = "h"),
       filename = paste0(figurePath, "Fig_4_FRSC.png"), 
       dpi = dpi,
       width = 90,
       height = 45,
       units = "mm")

# Save as .svg
ggsave(plot_grid(p1_mae_draws, pred_means_subject_draws_p, align = "h"),
       filename = paste0(figurePath, "Fig_4_FRSC.svg"), 
       dpi = dpi,
       width = 90,
       height = 45,
       units = "mm")

Report MAE

  • Average marginal effect -0.08 SDs (95 % CI [-0.11, -0.06])

Methods

Task and virtual environment

Trial types of the experiment

Encoding trial

Average trial lengths:

# Subset data to encoding only
encoding_only <- OLM_7T_trial_results[OLM_7T_trial_results$trialType == "encoding", ]

# Calculate trial duration
encoding_only$trial_duration <- encoding_only$end_time - encoding_only$start_time

# Calculate the average for each subject
trial_dur <- ddply(encoding_only, c("subject"), summarise,
                   total_length = sum(trial_duration),
                   run_length = total_length/2,
                   trial_duration = mean(trial_duration),
                   N = length(subject))

# Create report strings
str1   <-  mean_SD_str2(trial_dur$trial_duration, report_type, digits1, rounding_type)
str2   <-  mean_SD_str2(trial_dur$run_length, report_type, digits1, rounding_type)

Average trial duration of an encoding trial: M = 20 (SD = 3.1).

Average encoding run duration: M = 360 (SD = 55).

Procedure

Days between sessions

# Load lab notebook with information on each subject
lab_notebook <- read_excel("data/ignore_fMRI_version1/VP00157_PPT_Info_Final.xlsx")

# Subset to only the participants used in this analysis
lab_notebook <- lab_notebook[lab_notebook$`ID Number` %in% subjects, ]

# Extract time between sessions
val <- lab_notebook$`Delay between sessions`

# Create report string
str1   <-  mean_SD_str2(val, report_type, digits1, rounding_type)
  • Number of days between sessions M = 12 (SD = 13)

MRI acquisition

Number of volumes

# Load the .RData with the values
load("data/ignore_fMRI_version1/extracted_values/SpaNov_nVol_mean.RData")

# Calculate the total number of volumes
nVol_mean$nVol_total <- nVol_mean$nVol1 + nVol_mean$nVol2

# Convert this to minutes
nVol_mean$dataMinutes <- nVol_mean$nVol_total/60

# Create report strings
str1   <-  mean_SD_str2(nVol_mean$nVol1, report_type, digits1, rounding_type)
str2   <-  mean_SD_str2(nVol_mean$nVol2, report_type, digits1, rounding_type)
minTime <- signif(min(nVol_mean$dataMinutes), 3)
  • Number of volumes Run 1: M = 410 (SD = 81)
  • Number of volumes Run 2: M = 370 (SD = 37)
  • Minimum data in minutes: 10.5

Neuroimaging analysis overview

Standard GLMs

Gradient analysis

Hippocampal gradient analysis

Get the number of average voxels per position

voxelNum <- mean_SD_str2(HC_data_agg_pos$n, 1, digits1, rounding_type)

The average number of voxels per position is M = 7.6 (SD = 3)

Supplementary material

Supplementary notes

Behavioural models of novelty separately

# Visits
load("fitted_brms_models/SpaNov_m_visits_poisson.Rdata")
mae1 <- avg_comparisons(m_visits_run1, variables = "s_onset_rel", type = "response", re_formula = NULL)
mae2 <- avg_comparisons(m_visits_run2, variables = "s_onset_rel", type = "response", re_formula = NULL)

# Time since
load("fitted_brms_models/SpaNov_m_timeSince.Rdata")
mae3 <- avg_comparisons(m_timeSince_run1, variables = "s_onset_rel", type = "response", re_formula = NULL)
mae4 <- avg_comparisons(m_timeSince_run2, variables = "s_onset_rel", type = "response", re_formula = NULL)

# Make strings
str1 <- create_str_from_avg_comparisons(mae1, "visits")
str2 <- create_str_from_avg_comparisons(mae2, "visits")
str3 <- create_str_from_avg_comparisons(mae3, "seconds")
str4 <- create_str_from_avg_comparisons(mae4, "seconds")

Average marginal effects: - Run 1: 1.85 visits (95 % CI [1.63, 2.09]) - Run 2: 2.16 visits (95 % CI [1.95, 2.37]) - Run 1: 56.34 seconds (95 % CI [49.46, 63.84]) - Run 2: -145.42 seconds (95 % CI [-157.7, -133.73])

Supplementary tables

Significant clusters for linear contrast (Level 1 < Level 2 < Level 3 < …)

Click here for table.
part cluster size zValue_mean zValue_median zValue_max zValue_min zValue_peak_abs cluster_mass clsuter_mass_log peak_x peak_y peak_z peak_region peak_region_name num_regions
cortex_left 1 22 3.063706 2.999546 3.6 2.5 3.638579 67 4.210668 -3.791392 -46.5213089 30.8302307 214 d23ab 2
cortex_left 2 479 3.210936 3.128089 4.6 2.5 4.618636 1500 7.338263 -22.350677 -43.2125931 72.1212463 232 2 7
cortex_left 3 357 3.129615 3.080705 4.4 2.5 4.356853 1100 7.018646 -44.840401 -25.9435463 17.5734921 284 RI 8
cortex_left 4 1425 3.630920 3.513689 6.4 2.5 6.428157 5200 8.551413 -5.515844 -16.9644909 58.8415260 220 24dd 16
cortex_left 9 30 2.695159 2.697803 2.9 2.5 2.874365 81 4.392655 -9.651464 -42.0632935 60.6421967 217 5mv 2
cortex_left 12 468 3.346280 3.275896 5.1 2.5 5.148257 1600 7.356318 -51.256317 -17.4799042 51.5853081 189 3b 6
cortex_left 18 44 2.766905 2.739925 3.1 2.5 3.093434 120 4.801919 -59.008034 -37.4201393 -0.5577022 309 STSdp 2
cortex_left 19 70 3.043194 2.963094 3.9 2.6 3.897602 210 5.361403 -61.696888 -43.7761002 11.0784178 319 TPOJ1 2
cortex_left 20 28 2.984275 2.937808 3.5 2.5 3.458510 84 4.425561 -61.709724 -46.8277969 -2.5455508 313 TE1p 2
cortex_left 24 46 3.100437 3.043691 3.7 2.5 3.705008 140 4.960185 -44.560581 -14.6814651 0.5779716 283 52 4
cortex_left 26 272 3.135538 2.986855 4.8 2.5 4.801907 850 6.748603 -58.528763 9.8367786 16.2200336 258 6r 8
cortex_left 39 109 3.257279 3.251825 4.2 2.6 4.240385 360 5.872240 -50.125221 -73.3484497 8.6767569 182 MST 4
cortex_left 40 180 3.216062 3.199615 4.4 2.5 4.369700 580 6.361115 -51.717174 -66.0832825 29.1244259 330 PGi 2
cortex_left 42 23 2.880656 2.844342 3.3 2.5 3.317583 66 4.193512 -54.612988 -0.3613958 42.1993637 236 6v 4
cortex_left 43 77 2.885721 2.847735 3.4 2.5 3.436671 220 5.403580 -62.776825 -33.4391670 32.7153397 328 PF 3
cortex_left 47 121 2.972278 2.908775 3.8 2.4 3.818675 360 5.885119 -7.869712 47.3418388 -12.9690065 245 10r 4
cortex_left 55 150 3.014078 2.915844 4.0 2.5 4.010159 450 6.113929 -55.074410 4.9239883 -23.9912205 308 STSda 5
cortex_left 65 279 3.013078 2.995985 3.8 2.4 3.759457 840 6.734174 -13.645574 46.4050636 47.3006058 251 9p 7
cortex_left 68 29 3.141467 3.073398 4.0 2.5 4.037845 91 4.511986 -8.394770 51.6081810 7.1561594 249 9m 3
cortex_left 71 42 2.948894 2.914134 3.6 2.5 3.602218 120 4.819100 -34.522411 49.6495667 31.6823597 264 46 1
cortex_right 87 1643 3.457321 3.398987 5.7 2.4 5.702450 5700 8.644773 13.677155 -16.6584091 75.3533707 55 6mp 20
cortex_right 92 505 3.332931 3.249488 5.0 2.5 5.003071 1700 7.428411 33.398880 -40.5589371 66.4627914 52 2 6
cortex_right 93 31 2.747490 2.681356 3.2 2.5 3.156602 85 4.444675 30.130425 -26.7811832 67.4318771 53 3a 2
cortex_right 95 20 2.676820 2.652724 2.9 2.5 2.918891 54 3.980362 22.743727 -42.0650024 73.2298584 52 2 1
cortex_right 96 141 3.187271 3.139549 4.2 2.5 4.176615 450 6.107925 57.448917 -13.8404722 48.2410164 51 1 3
cortex_right 99 457 3.089366 2.971921 4.8 2.4 4.844986 1400 7.252649 47.416145 -25.7598839 18.3474407 104 RI 9
cortex_right 102 80 3.276918 3.268695 4.5 2.6 4.506863 260 5.568930 62.873524 -40.0575027 15.5642920 28 STV 3
cortex_right 103 37 2.822053 2.749717 3.4 2.4 3.436898 100 4.648383 57.842972 -34.0697250 -1.2721851 129 STSdp 2
cortex_right 110 21 2.753690 2.766122 3.2 2.5 3.189069 58 4.057464 46.799515 -2.4484978 -9.9124584 167 PoI1 2
cortex_right 112 55 3.004002 2.978083 3.7 2.6 3.667743 170 5.107278 45.115521 8.0075855 8.4428720 115 FOP2 2
cortex_right 113 25 2.951910 2.927984 3.4 2.5 3.383995 74 4.301328 41.860241 29.4441929 3.9625542 108 FOP4 2
cortex_right 116 37 3.206650 3.216331 3.9 2.4 3.909213 120 4.776145 16.628929 -83.4581909 40.5141487 3 V6 2
cortex_right 117 220 3.506168 3.399297 5.7 2.4 5.696318 770 6.648151 50.564453 -72.4114990 6.8142180 23 MT 7
cortex_right 125 216 3.174723 3.072107 4.3 2.5 4.313631 690 6.530499 61.787083 8.8548021 19.2581615 56 6v 6
cortex_right 130 24 2.946266 2.877828 3.7 2.5 3.672417 71 4.258592 43.409565 17.7264061 -35.5937271 131 TGd 1
cortex_right 146 38 2.828438 2.834183 3.3 2.5 3.254191 110 4.677311 8.638423 55.7692184 25.6142464 69 9m 2
cortex_right 147 25 2.989004 2.958746 3.6 2.5 3.600116 75 4.313816 35.940197 50.5307732 30.8277359 84 46 1
cortex_right 155 26 2.818180 2.809673 3.1 2.6 3.136499 73 4.294188 58.314800 0.0648438 -16.3995800 128 STSda 1
subcort 168 161 3.074034 2.976389 4.5 2.4 4.542136 490 6.204395 26.000000 -88.0000000 -38.0000000 371 CEREBELLUM_RIGHT 1
subcort 193 107 2.998622 2.911509 4.1 2.5 4.123521 320 5.770982 20.000000 -14.0000000 -20.0000000 376 HIPPOCAMPUS_RIGHT 2
subcort 194 135 3.144796 3.095737 4.5 2.5 4.509548 420 6.051024 -18.000000 -12.0000000 -20.0000000 367 HIPPOCAMPUS_LEFT 2
subcort 215 43 2.861746 2.787239 3.8 2.4 3.830762 120 4.812632 -16.000000 8.0000000 -8.0000000 364 PUTAMEN_LEFT 1
part cluster regions labels
cortex_left 1 215,214 31pv,d23ab
cortex_left 2 232,231,227,222,229,219,225 2,1,7PC,7AL,VIP,5L,7Am
cortex_left 3 282,328,285,284,204,205,354,281 OP2-3,PF,PFcm,RI,A1,PSL,LBelt,OP1
cortex_left 4 218,221,223,240,237,239,189,233,188,235,224,216,217,220,234,276 23c,24dv,SCEF,p32pr,p24pr,a24pr,3b,3a,4,6mp,6ma,5m,5mv,24dd,6d,6a
cortex_left 9 219,217 5L,5mv
cortex_left 12 189,233,188,190,234,231 3b,3a,4,FEF,6d,1
cortex_left 18 310,309 STSvp,STSdp
cortex_left 19 319,208 TPOJ1,STV
cortex_left 20 310,313 STSvp,TE1p
cortex_left 24 282,283,347,348 OP2-3,52,PoI1,Ig
cortex_left 26 288,295,294,293,258,279,236,188 FOP4,FOP2,FOP3,FOP1,6r,43,6v,4
cortex_left 39 203,182,320,321 MT,MST,TPOJ2,TPOJ3
cortex_left 40 330,331 PGi,PGs
cortex_left 42 188,236,191,192 4,6v,PEF,55b
cortex_left 43 328,327,285 PF,PFop,PFcm
cortex_left 47 268,245,345,252 10v,10r,s32,10d
cortex_left 55 311,303,308,356,312 TGd,STGa,STSda,STSva,TE1a
cortex_left 65 206,250,251,252,249,267,248 SFL,8BL,9p,10d,9m,9a,8Ad
cortex_left 68 241,244,249 a24,p32,9m
cortex_left 71 264 46
cortex_right 87 12,38,41,43,60,57,8,55,44,37,39,36,40,52,51,10,9,54,53,96 55b,23c,24dv,SCEF,p32pr,p24pr,4,6mp,6ma,5mv,5ROI,5m,24dd,2,1,FEF,3b,6d,3a,6a
cortex_right 92 39,52,47,42,49,48 5ROI,2,7PC,7AROI,VIP,LIPv
cortex_right 93 9,53 3b,3a
cortex_right 95 52 2
cortex_right 96 51,52,116 1,2,PFt
cortex_right 99 148,105,104,25,124,174,101,147,100 PF,PFcm,RI,PSROI,PBelt,LBelt,OP1,PFop,OP4
cortex_right 102 139,28,25 TPOJ1,STV,PSROI
cortex_right 103 129,125 STSdp,A5
cortex_right 110 167,178 PoI1,PI
cortex_right 112 115,114 FOP2,FOP3
cortex_right 113 108,169 FOP4,FOP5
cortex_right 116 152,3 V6A,V6
cortex_right 117 156,2,23,137,157,140,141 V4t,MST,MT,PHT,FST,TPOJ2,TPOJ3
cortex_right 125 113,78,99,56,11,12 FOP1,6r,43,6v,PEF,55b
cortex_right 130 131 TGd
cortex_right 146 87,69 9a,9m
cortex_right 147 84 46
cortex_right 155 128 STSda
subcort 168 371 CEREBELLUM_RIGHT
subcort 193 376,377 HIPPOCAMPUS_RIGHT,AMYGDALA_RIGHT
subcort 194 367,368 HIPPOCAMPUS_LEFT,AMYGDALA_LEFT
subcort 215 364 PUTAMEN_LEFT

Significant clusters for linear contrast (Level 1 > Level 2 > Level 3 > …)

Click here for table.
part cluster size zValue_mean zValue_median zValue_max zValue_min zValue_peak_abs cluster_mass clsuter_mass_log peak_x peak_y peak_z peak_region peak_region_name num_regions
cortex_left 1 90 -3.519341 -3.443073 -2.4 -4.8 4.790668 320 5.758084 -2.259454 -31.05509 32.052841 194 RSC 3
cortex_left 2 2060 -3.930611 -3.876594 -2.5 -6.7 6.735095 8100 8.999256 -32.472549 -88.98758 29.779585 326 IP0 28
cortex_left 3 221 -3.194249 -3.094511 -2.5 -4.3 4.344084 710 6.559515 -28.574898 -52.22017 -18.019115 307 PHA3 10
cortex_left 5 27 -2.856392 -2.848945 -2.5 -3.3 3.321889 77 4.345396 -8.653727 -54.54932 5.236495 211 POS1 2
cortex_left 7 84 -3.043542 -3.026831 -2.5 -3.7 3.736013 260 5.543839 -5.988591 25.43487 52.024731 243 8BM 4
cortex_left 11 216 -3.485387 -3.370426 -2.5 -5.1 5.099926 750 6.623857 -27.186953 19.61117 58.827522 277 i6-8 8
cortex_left 14 107 -3.212653 -3.204306 -2.5 -4.1 4.055283 340 5.839926 -47.295277 16.06033 34.457317 260 IFJp 4
cortex_left 15 131 -3.085281 -3.057481 -2.6 -3.9 3.863326 400 6.001840 -44.775909 41.61255 25.340578 263 p9-46v 4
cortex_left 17 55 -3.049246 -2.937280 -2.5 -3.8 3.841370 170 5.122228 -36.455795 60.54835 7.798744 265 a9-46v 3
cortex_left 29 483 -4.251315 -4.164235 -2.5 -6.8 6.760676 2100 7.627245 -13.782526 -92.53966 -12.004962 184 V2 4
cortex_left 32 22 -2.944452 -2.925727 -2.5 -3.4 3.412598 65 4.170965 -6.316596 -90.02340 8.658187 181 V1 1
cortex_right 38 2303 -3.929091 -3.748647 -2.4 -6.9 6.878045 9000 9.110376 36.958008 -87.18770 24.761272 146 IP0 35
cortex_right 39 244 -3.485530 -3.477159 -2.5 -4.9 4.867373 850 6.745788 24.893423 -51.38643 -16.057383 160 VMV2 7
cortex_right 40 118 -3.852320 -3.901859 -2.5 -4.9 4.903612 450 6.119360 3.001509 -31.90833 32.272301 14 RSC 1
cortex_right 44 63 -3.120948 -3.171131 -2.5 -3.8 3.839579 200 5.281271 6.831443 29.15144 44.001465 179 a32pr 3
cortex_right 46 402 -3.335551 -3.302893 -2.5 -4.7 4.685742 1300 7.201090 31.845680 21.57755 58.013378 97 i6-8 10
cortex_right 50 186 -3.224232 -3.150901 -2.5 -4.5 4.506076 600 6.396442 48.756924 39.91827 20.827803 81 IFSp 4
cortex_right 70 472 -4.438338 -4.449830 -2.4 -7.2 7.179883 2100 7.647259 13.675855 -91.91169 -6.089150 1 V1 5
cortex_right 74 28 -2.969049 -2.899777 -2.6 -3.6 3.633223 83 4.420446 10.256091 -91.39706 9.405240 1 V1 2
cortex_right 77 58 -3.031525 -2.926241 -2.5 -4.5 4.476549 180 5.169509 23.169453 66.53020 2.306681 170 p10p 3
cortex_right 80 23 -2.992635 -3.061711 -2.5 -3.5 3.546734 69 4.231649 26.927168 42.81836 44.168594 68 8Ad 2
subcort 82 172 -3.229611 -3.108652 -2.5 -4.8 4.768074 560 6.319856 -38.000000 -70.00000 -52.000000 361 CEREBELLUM_LEFT 1
subcort 84 283 -3.336223 -3.222986 -2.4 -5.3 5.324766 940 6.850286 10.000000 -46.00000 -48.000000 371 CEREBELLUM_RIGHT 2
subcort 86 185 -3.337890 -3.209643 -2.4 -4.7 4.740293 620 6.425695 34.000000 -70.00000 -50.000000 371 CEREBELLUM_RIGHT 1
subcort 87 198 -3.299109 -3.250997 -2.5 -4.8 4.785460 650 6.481920 -24.000000 -32.00000 -42.000000 361 CEREBELLUM_LEFT 1
subcort 92 372 -3.564598 -3.446843 -2.5 -5.7 5.698634 1300 7.189945 -4.000000 -80.00000 -40.000000 361 CEREBELLUM_LEFT 1
subcort 99 271 -3.506712 -3.379627 -2.5 -6.0 6.013061 950 6.856798 8.000000 -78.00000 -24.000000 371 CEREBELLUM_RIGHT 1
subcort 106 20 -3.091158 -3.096199 -2.6 -4.0 4.030955 62 4.124278 -10.000000 -58.00000 -38.000000 361 CEREBELLUM_LEFT 1
subcort 112 201 -3.193322 -3.150749 -2.5 -4.6 4.587187 640 6.464367 28.000000 -64.00000 -30.000000 371 CEREBELLUM_RIGHT 1
subcort 113 250 -3.332495 -3.219638 -2.5 -5.0 4.964678 830 6.725182 -26.000000 -60.00000 -32.000000 361 CEREBELLUM_LEFT 1
subcort 141 47 -2.979207 -2.928996 -2.5 -4.0 4.010671 140 4.941805 18.000000 -32.00000 6.000000 372 THALAMUS_RIGHT 1
part cluster regions labels
cortex_left 1 194,214,212 RSC,d23ab,23d
cortex_left 2 211,301,207,218,342,297,324,329,230,197,195,322,226,225,209,210,199,326,323,331,325,275,228,186,338,321,330,184 POS1,ProS,PCV,23c,31a,AIP,IP2,PFm,MIP,IPS1,POS2,DVT,7PL,7Am,7Pm,7m,V3B,IP0,PGp,PGs,IP1,LIPd,LIPv,V4,V3CD,TPOJ3,PGi,V2
cortex_left 3 306,333,301,299,184,340,307,335,343,334 PHA1,VMV1,ProS,PreS,V2,VMV2,PHA3,PHA2,VVC,VMV3
cortex_left 5 194,211 RSC,POS1
cortex_left 7 206,243,223,359 SFL,8BM,SCEF,a32pr
cortex_left 11 190,276,224,278,248,277,247,192 FEF,6a,6ma,s6-8,8Ad,i6-8,8Av,55b
cortex_left 14 253,259,260,191 8C,IFJa,IFJp,PEF
cortex_left 15 253,263,262,261 8C,p9-46v,IFSa,IFSp
cortex_left 17 265,266,350 a9-46v,9-46d,p10p
cortex_left 29 181,184,185,186 V1,V2,V3,V4
cortex_left 32 181 V1
cortex_right 38 149,31,121,14,162,27,35,38,34,144,117,50,17,30,15,142,46,45,152,29,161,6,19,146,143,151,13,145,95,48,20,158,150,1,3 PFm,POS1,ProS,RSC,31a,PCV,31pv,23c,d23ab,IP2,AIP,MIP,IPS1,7m,POS2,DVT,7PROI,7Am,V6A,7Pm,31pd,V4,V3B,IP0,PGp,PGs,V3A,IP1,LIPd,LIPv,LO1,V3CD,PGi,V1,V6
cortex_right 39 160,126,153,127,155,163,154 VMV2,PHA1,VMV1,PHA3,PHA2,VVC,VMV3
cortex_right 40 14 RSC
cortex_right 44 63,43,179 8BM,SCEF,a32pr
cortex_right 46 10,96,73,80,79,11,98,68,97,67 FEF,6a,8C,IFJp,IFJa,PEF,s6-8,8Ad,i6-8,8Av
cortex_right 50 73,83,81,82 8C,p9-46v,IFSp,IFSa
cortex_right 70 1,4,5,6,7 V1,V2,V3,V4,V8
cortex_right 74 1,4 V1,V2
cortex_right 77 170,90,72 p10p,10pp,10d
cortex_right 80 68,84 8Ad,46
subcort 82 361 CEREBELLUM_LEFT
subcort 84 371,366 CEREBELLUM_RIGHT,BRAIN_STEM
subcort 86 371 CEREBELLUM_RIGHT
subcort 87 361 CEREBELLUM_LEFT
subcort 92 361 CEREBELLUM_LEFT
subcort 99 371 CEREBELLUM_RIGHT
subcort 106 361 CEREBELLUM_LEFT
subcort 112 371 CEREBELLUM_RIGHT
subcort 113 361 CEREBELLUM_LEFT
subcort 141 372 THALAMUS_RIGHT

Supplementary figures

Learning curves

# Creating these figures takes too long so eval = FALSE
load("fitted_brms_models/SpaNov_m_navTime.Rdata")

# Subset to retrieval trials
OLM_7T_encoding <- OLM_7T_trial_results[OLM_7T_trial_results$trialType == "encoding", ]

# Calculate distance between start and object
x1 <- OLM_7T_encoding$start_x
z1 <- OLM_7T_encoding$start_z
x2 <- OLM_7T_encoding$object_x
z2 <- OLM_7T_encoding$object_z
OLM_7T_encoding$start2obj <- euclideanDistance3D(x1, 1, z1, x2, 1, z2)

# Scale values
OLM_7T_encoding$s_timesObjectPresented <- scale_m0_sd0.5(OLM_7T_encoding$timesObjectPresented)
OLM_7T_encoding$s_start2obj <- scale_m0_sd0.5(OLM_7T_encoding$start2obj)

# Rename model & set variable of interest
m   <- m_navTime3
voi <- "s_timesObjectPresented"

# Get which raw values the scaled values correspond to
raw_and_scaled_values <- ddply(OLM_7T_encoding, c("timesObjectPresented"), summarise,
                          raw = timesObjectPresented[1],
                          scaled = mean(s_timesObjectPresented))

# Calculate marginal average effect is the contrast 0 vs. 1
mae1_navTime_respone <- avg_comparisons(m, variables = voi, type = "response", re_formula = NULL)
#mae1_navTime_link    <- avg_comparisons(m, variables = voi, type = "link", re_formula = NULL)

# Get posterior draws
mae1_navTime_respone_draws <-  get_draws(mae1_navTime_respone)

p1_mae_draws <- ggplot(mae1_navTime_respone_draws, aes(x = draw)) +
  geom_vline(xintercept = 0, linetype = "dashed") + 
  stat_halfeye(fill = base_col[1], point_size = 1.5) +
  labs(x = "Navigation time (seconds)", y = "Density", 
       title = "Marginal average effect\n of number of times presented") +
  theme_journal() +
  theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5)) 

# Create estimates for line
grid <- datagrid(s_timesObjectPresented = seq_min_2_max(raw_and_scaled_values$scaled), 
                 subject = unique, s_start2obj = 0, model = m)
pred_line_navTime <- predictions(m, newdata = grid, by = voi, re_formula = NULL, type = "response")
pred_line_navTime <- as.data.frame(pred_line_navTime)

# Create estimates actually existing values
grid <- datagrid(s_timesObjectPresented = unique, subject = unique, s_start2obj = 0, model = m)
pred_actual_navTime <- predictions(m, newdata = grid, by = voi, re_formula = NULL, type = "response")
pred_actual_navTime <- as.data.frame(pred_actual_navTime) 

# Calculate subject-level lines
grid <- datagrid(s_timesObjectPresented = seq_min_2_max(raw_and_scaled_values$scaled), 
                 subject = unique, s_start2obj = 0, model = m)
pred_line_subject_navTime <- predictions(m, newdata = grid, type = "response")
pred_line_subject_navTime <- as.data.frame(pred_line_subject_navTime)

# Make DFs unique
pred_line_navTime1         <- pred_line_navTime
pred_actual_navTime1       <- pred_actual_navTime
pred_line_subject_navTime1 <- pred_line_subject_navTime

# Assemble plot
p1_line <-ggplot() +
  geom_line(data = pred_line_subject_navTime1, mapping = aes(x = s_timesObjectPresented, y = estimate, group = subject),
            colour = "darkgrey") +
  geom_line(data = pred_line_navTime1, mapping = aes(x = s_timesObjectPresented, y = estimate), linewidth = 0.3) +
  geom_point(data = pred_actual_navTime1, mapping = aes(x = s_timesObjectPresented, y = estimate)) +
  geom_errorbar(data = pred_actual_navTime1, mapping = aes(x = s_timesObjectPresented, ymin = conf.low, ymax = conf.high),
                 width = 0.02) +
  scale_x_continuous(breaks = raw_and_scaled_values$scaled, labels = raw_and_scaled_values$raw) +
  theme_journal() +
  labs(x = "Number of times the object was presented", y = "Navigation time\n(seconds)")

###########################################################################_
# Assemble figure
learning_curves <- plot_grid(p1_mae_draws, p1_line,  align = "hv", ncol = 2, 
                             labels = "auto", label_size = 10)

# Save 
ggsave(learning_curves,
       filename = paste0(figurePath, "Fig_S01_learning_curves.png"), 
       dpi = dpi, width = 130, height = 40, units = "mm")

Report MAEs

# Load the model 
load("fitted_brms_models/SpaNov_m_navTime.Rdata")

# Rename model & set variable of interest
m   <- m_navTime3
voi <- "s_timesObjectPresented"
mae1_navTime_respone <- avg_comparisons(m, variables = voi, type = "response", re_formula = NULL)
str1 <- create_str_from_avg_comparisons(mae1_navTime_respone, "seconds")
  • Average marginal effect: str 1

Locomotion by novelty

# Load the image from the gradient event file creation
load("event_tables/images/SpaNov_event_file_gradients.RData")

# Create a data frame
ev_info <- rbindlist(condition_info)

# Subset to the analysis with 6 levels and encoding and for the novelty score
ev_info <- ev_info[ev_info$curr_numLvl == 6 & ev_info$runType == "e" & ev_info$type == "noveltyScore", ]

# Calculate average across runs
ev_info_agg <- ddply(ev_info, c("subj", "level"), summarise, 
                     mean_per_tra_arc = mean(mean_per_tra_arc),
                     mean_per_rot_arc = mean(mean_per_rot_arc),
                     mean_per_sta_arc = mean(mean_per_sta_arc))

# Calculate average for each run
ev_info_agg1 <- ddply(ev_info, c("subj", "run"), summarise, 
                     mean_per_tra_arc = mean(mean_per_tra_arc),
                     mean_per_rot_arc = mean(mean_per_rot_arc),
                     mean_per_sta_arc = mean(mean_per_sta_arc))


# Remove middle levels
ev_info_agg <- ev_info_agg[ev_info_agg$level == "lvl1" | ev_info_agg$level == "lvl6", ]

# Create the plots
## Translation
p1 <- ggplot(ev_info_agg, aes(x = level, y = mean_per_tra_arc, fill = level)) +
  geom_line(aes(group = subj), colour = boxplot_line_col) +
  geom_point(colour = boxplot_point_col) +
  geom_boxplot(colour = boxplot_border_col, outlier.shape = NA, width = 0.4) +
  stat_summary(geom = "point", fun = "mean", col = 'white', shape = 24, 
               position = position_dodge(width =  0.75), size = 0.8) +
  scale_fill_manual(values = novFam_gradient[c(1, 6)]) +
  theme_journal() + 
  scale_y_continuous(breaks = c(0.25, 0.82, 1.40)) + 
  coord_cartesian(xlim = c(0.5, 2.5), ylim = c(0.25, 1.40), expand = FALSE) +
  theme(legend.position = "none") +
  labs(title = "Translation", x = "Novelty level", y = "arcsine(percentage)")

## Rotation
p2 <- ggplot(ev_info_agg, aes(x = level, y = mean_per_rot_arc, fill = level)) +
  geom_line(aes(group = subj), colour = boxplot_line_col) +
  geom_point(colour = boxplot_point_col) +
  geom_boxplot(colour = boxplot_border_col, outlier.shape = NA, width = 0.4) +
  stat_summary(geom = "point", fun = "mean", col = 'white', shape = 24, 
               position = position_dodge(width =  0.75), size = 0.8) +
  scale_fill_manual(values = novFam_gradient[c(1, 6)]) +
  theme_journal() + 
  scale_y_continuous(breaks = c(-1.5, -1, -0.58)) + 
  coord_cartesian(xlim = c(0.5, 2.5), ylim = c(-1.5, -0.58), expand = FALSE) +
  theme(legend.position = "none") +
  labs(title = "Rotation", x = "Novelty level", y = "arcsine(percentage)")

## Stationary
p3 <- ggplot(ev_info_agg, aes(x = level, y = mean_per_sta_arc, fill = level)) +
  geom_line(aes(group = subj), colour = boxplot_line_col) +
  geom_point(colour = boxplot_point_col) +
  geom_boxplot(colour = boxplot_border_col, outlier.shape = NA, width = 0.4) +
  stat_summary(geom = "point", fun = "mean", col = 'white', shape = 24, 
               position = position_dodge(width =  0.75), size = 0.8) +
  scale_fill_manual(values = novFam_gradient[c(1, 6)]) +
  theme_journal() + 
  scale_y_continuous(breaks = c(-1.5, -1.10, -0.73)) + 
  coord_cartesian(xlim = c(0.5, 2.5), ylim = c(-1.5, -0.73), expand = FALSE) +
  theme(legend.position = "none") +
  labs(title = "Stationary", x = "Novelty level", y = "arcsine(percentage)")

# Combine and save
p_comb <- plot_grid(p1, p2, p3, ncol = 3, align = "hv")

ggsave(p_comb,
       filename = paste0(figurePath, "Fig_S02_locomotion.png"), 
       dpi = dpi,
       width = 90,
       height = 40,
       units = "mm")

# Calculate the test statistics
## Translation
val1     <- ev_info_agg[ev_info_agg$level == 'lvl1', "mean_per_tra_arc"]
val2     <- ev_info_agg[ev_info_agg$level == 'lvl6', "mean_per_tra_arc"]
diff_val <- val1 - val2
BF1      <- signif(reportBF(ttestBF(diff_val)), digits1)
d1       <- signif(mean(diff_val)/sd(diff_val), digits1)

## Rotation
val1     <- ev_info_agg[ev_info_agg$level == 'lvl1', "mean_per_rot_arc"]
val2     <- ev_info_agg[ev_info_agg$level == 'lvl6', "mean_per_rot_arc"]
diff_val <- val1 - val2
BF2      <- signif(reportBF(ttestBF(diff_val)), digits1)
d2       <- signif(mean(diff_val)/sd(diff_val), digits1)

## Stationary
val1     <- ev_info_agg[ev_info_agg$level == 'lvl1', "mean_per_sta_arc"]
val2     <- ev_info_agg[ev_info_agg$level == 'lvl6', "mean_per_sta_arc"]
diff_val <- val1 - val2
BF3      <- signif(reportBF(ttestBF(diff_val)), digits1)
d3       <- signif(mean(diff_val)/sd(diff_val), digits1)
  • Translation Level 1 vs. Level 6: BF10 = 3.5^{15}, d = -1.7
  • Rotation Level 1 vs. Level 6: BF10 = 7.3^{16}, d = 1.9
  • Stationary Level 1 vs. Level 6: BF10 = 9.4^{8}, d = 1.1

Spatial distribution of novelty scores

# Circle definitions
## Area of the arena
arenaDiameter <- 180
radius1       <- arenaDiameter/2
areaCircle1   <- pi*(radius1^2)

## Area of circle with half the radius
radius2     <- radius1/2
areaCircle2 <- pi*radius2^2

# Draw inner and outer circle
## Outer circle
circle1   <- circleFun(c(0, 0), arenaDiameter, npoints = 1000)
circle1$z <- circle1$y # Rename to avoid error that columns do not match
circle1$y <- NULL

## Inner circle
circle2   <- circleFun(c(0, 0), radius2*2, npoints = 1000)
circle2$z <- circle2$y # Rename to avoid error that columns do not match
circle2$y <- NULL

# Convert list to data frame
event_info_df <- rbindlist(event_info, idcol = "id")
event_info_df <- as.data.frame(event_info_df)

# Subset to gradient level 1
event_info_df_lvl1 <- event_info_df[event_info_df$noveltyScore_gradientLevel == "lvl1" & 
                                    !is.na(event_info_df$noveltyScore_gradientLevel), ]
event_info_df_lvl1 <- event_info_df_lvl1[event_info_df_lvl1$run %in% c(1, 3), ]

spa_nov_dist <- ggplot() +
  # Circles
  geom_polygon(data = circle1, mapping = aes(x = x, y = z), fill = "lightgrey") +
  geom_polygon(data = circle2, mapping = aes(x = x, y = z), fill = "darkgrey") +
  # Lines and text for outer radius
  geom_segment(aes(x = 0, y = -97, xend = radius1, yend = -97)) +
  geom_segment(aes(x = radius1, y = 0, xend = radius1, yend = -97), linetype = "dashed") +
  geom_segment(aes(x = 0, y = 0, xend = 0, yend = -97), linetype = "dashed") +
  annotate("text", x = 45, y = -102, label = c("90 vm"), size = 2) +
  # Lines for inner radius
  geom_segment(aes(x = 0, y = -95, xend = radius2, yend = -95)) +
  geom_segment(aes(x = radius2, y = 0, xend = radius2, yend = -95), linetype = "dashed") +
  annotate("text", x = radius2/2, y = -90, label = c("45 vm"), size = 2) +
  # Points
  geom_point(data = event_info_df_lvl1, mapping = aes(x = pos_x, y = pos_z)) +
  # Arrows and text
  geom_curve(aes(x = -95, y = 70, xend = -50, yend = 50), ncp = 10,  arrow = arrow(length = unit(1, "mm"))) +
  geom_curve(aes(x = 95, y = 70, xend = 30, yend = 0), ncp = 10,  arrow = arrow(length = unit(1, "mm")), curvature = -0.5) +
  annotate("text", x = c(-95, 95), y = c(73, 73), label = c("Periphery", "Centre"), size = 2) +
  # Other settings
  coord_equal(xlim = c(-105, 105)) +
  theme_void() +
  theme(plot.background = element_rect(fill = "white", colour = "white"))

ggsave(spa_nov_dist,
       filename = paste0(figurePath, "Fig_S03_spa_nov_dist.png"), 
       dpi = dpi, width = 80, height = 80, units = "mm")

# Calculate centrality for each level
## Subset to encoding runs & none NA values
event_info_df_encoding <- event_info_df[event_info_df$run %in% c(1, 3), ]
event_info_df_encoding <- na.omit(event_info_df_encoding)

## Calculate distance to centre
event_info_df_encoding$dist2centre <- euclideanDistance3D(0, 1, 0, event_info_df_encoding$pos_x, 1, event_info_df_encoding$pos_z)

## Categorise centrality
event_info_df_encoding$centrality1 <- ifelse(event_info_df_encoding$dist2centre  >= radius2, 'perihery', 'central')

## Calculate average for each level
avg_level_centrality <- ddply(event_info_df_encoding, c("noveltyScore_gradientLevel"), summarise,
                              centrality1_score = mean(centrality1 == "central"))

Minimal influence of centrality on spatial novelty effects

# Loading files
dist2centre_file <- "/media/alex/work/Seafile/imaging_results/SpaNov/OLMe_7T_SpaNov_gradient_dist2centre-corrected_6lvl_smo4_MSMAll/cope7.feat/stats/vwc/results_lvl2cope1_dat_ztstat_c1.dscalar.nii"
SpaNov1_file     <- "/media/alex/work/Seafile/imaging_results/SpaNov/OLMe_7T_SpaNov_gradient_6lvl_cue-delay_smo4_MSMAll/cope7.feat/stats/vwc/results_lvl2cope1_dat_ztstat_c1.dscalar.nii"
dist2centre_xii  <- read_cifti(dist2centre_file)
SpaNov1_xii      <- read_cifti(SpaNov1_file)

# 
df <- data.frame(dist2centre = as.matrix(dist2centre_xii),
                 SpaNov1 = as.matrix(SpaNov1_xii))

# Create correlation string
r_val <- round(cor.test(df$dist2centre, df$SpaNov1)$estimate, 2)
p_val <- pValue(cor.test(df$dist2centre, df$SpaNov1)$p.value)
report_str <- paste0("r = ", r_val, ", p ", p_val)

# Create scatter plot
p_map_corr <- ggplot(df, aes(x = dist2centre, y = SpaNov1)) + 
  geom_density_2d_filled() +
  #geom_smooth(method = "lm", formula = "y ~ x", colour = "white", se = FALSE) + 
  theme_journal() +
  theme(legend.position = "none") +
  theme(plot.margin = unit(c(1, 4, 1, 1), "mm")) +
  coord_cartesian(expand = FALSE) +
  annotate("text", x = I(0.3), y = I(0.9), label = report_str, size = 2, colour = "white") +
  labs(x = "Z-map corrected for centrality", y = "Original z-map")

##############################################################################_
# Calculate centrality
## Area of the arena
radius1     <- 180/2
areaCircle1 <- pi*(radius1^2)

## Area of circle with half the radius
radius2     <- radius1/2
areaCircle2 <- pi*radius2^2

## Subset position data only encoding (include cue & delay because they are the same for all)
centrality_data <- locomotion_data

## Calculate distance to centre
centrality_data$dist2centre <- euclideanDistance3D(0, 1, 0, centrality_data$pos_x, 1, centrality_data$pos_z)

## Categorise centrality
centrality_data$centrality1 <- ifelse(centrality_data$dist2centre  >= radius2, 'perihery', 'central')

## Calculate average
avg_centrality <- ddply(centrality_data, c("ppid","subject"), summarise,
                        score1 = mean(centrality1 == "central"))

## Write .csv for PALM analysis
write.csv(avg_centrality[,-2], file = "data/ignore_fMRI_version1/SpaNov_centrality_designMatrix_data.csv",
          row.names = FALSE, quote = FALSE)

##############################################################################_
p <- plot_grid(p_map_corr, NULL, NULL, 
               nrow = 1, labels = "auto", label_size = 10, align = "hv")

# Combine plots
ggsave(p,
       filename = paste0(figurePath, "Fig_S04_control_for_centrality.png"), 
       dpi = dpi, width = 180, height = 57, units = "mm")

Individual differences in centrality do not predict spatial novelty

## Calculate average centrality
str1 <- mean_SD_str2(avg_centrality$score1, report_type, digits1, rounding_type)

## Create plot
avg_centrality_plot <- ggplot(avg_centrality, aes(x = score1)) +
  geom_histogram(colour = base_col[1], fill = base_col[1]) +  
  geom_vline(xintercept = mean(avg_centrality$score1), colour = "black", linetype = 2, size = 0.5) +
  scale_x_continuous(breaks = c(0.45, 0.55, 0.65)) + 
  scale_y_continuous(breaks = c(0, 5, 10)) + 
  coord_cartesian(xlim = c(0.44, 0.65), ylim = c(0, 10), 
                  expand = FALSE) +
  labs(x = "Centrality score", y = "Count") +
  theme_journal() +
  theme(plot.margin = unit(c(1, 2.5, 1, 2.5), "mm"))

##############################################################################_
#  Correlate with average beta values
## Add beta values to data frame
avg_centrality$WB_novelty     <- WB_novelty_values
avg_centrality$WB_familiarity <- WB_familiarity_values

## Calculate correlations
r1     <- round(cor(avg_centrality$score1, avg_centrality$WB_novelty), 2)
BF10_1 <- reportBF(correlationBF(avg_centrality$score1, avg_centrality$WB_novelty))
p1     <- pValue(cor.test(avg_centrality$score1, avg_centrality$WB_novelty)$p.value)
r2     <- round(cor(avg_centrality$score1, avg_centrality$WB_familiarity), 2)
BF10_2 <- reportBF(correlationBF(avg_centrality$score1, avg_centrality$WB_familiarity))
p2     <- pValue(cor.test(avg_centrality$score1, avg_centrality$WB_familiarity)$p.value)

scatter1 <- ggplot(avg_centrality, aes(x = score1, y = WB_novelty)) +
  geom_point() +
  geom_smooth(method = "lm", formula = "y ~ x", colour = base_col[1]) +
  #scale_x_continuous(breaks = c(0.50, 0.575, 0.65)) + 
  scale_y_continuous(breaks = c(-400, -200, 0, 200)) + 
  coord_cartesian(xlim = c(0.50, 0.65), ylim = c(-400, 200), 
                  expand = FALSE) +
  theme_journal() +
  labs(x = "Centrality score", y = "Whole-brain novelty\ncontrast (a.u.)")

scatter2 <- ggplot(avg_centrality, aes(x = score1, y = WB_familiarity)) +
  geom_point() +
  geom_smooth(method = "lm", formula = "y ~ x", colour = base_col[1]) +
  scale_y_continuous(breaks = c(-100, 0, 100, 200)) + 
  coord_cartesian(xlim = c(0.50, 0.65), ylim = c(-100, 210), 
                  expand = FALSE) +
  theme_journal() +
  labs(x = "Centrality score", y = "Whole-brain familiarity\ncontrast (a.u.)")

##############################################################################_
p <- plot_grid(avg_centrality_plot, scatter1, scatter2, 
               nrow = 1, labels = "auto", label_size = 10, align = "hv")

# Combine plots
ggsave(p,
       filename = paste0(figurePath, "Fig_S05_individual_differences_centrality.png"), 
       dpi = dpi, width = 180, height = 55, units = "mm")
  • The correlation between the centrality score and the novelty contrast was r = 0.07, BF10 = 0.35, p = .584.
  • The correlation between the centrality score and the familiarity contrast was r = 0.1, BF10 = 0.4, p = .442.
# Load the p-value maps that have been corrected across main effect and correlation
SpaNov1_centrality_p1_file <- "/media/alex/work/Seafile/imaging_results/SpaNov/OLMe_7T_SpaNov_gradient_6lvl_smo4_MSMAll_centrality/cope7.feat/stats/vwc/results_lvl2cope1_dat_ztstat_cfdrp_c3.dscalar.nii"
SpaNov1_centrality_p2_file <- "/media/alex/work/Seafile/imaging_results/SpaNov/OLMe_7T_SpaNov_gradient_6lvl_smo4_MSMAll_centrality/cope7.feat/stats/vwc/results_lvl2cope1_dat_ztstat_cfdrp_c4.dscalar.nii"

# Load the maps
SpaNov1_centrality_p1_xii <- read_cifti(SpaNov1_centrality_p1_file, brainstructures = "all")
SpaNov1_centrality_p2_xii <- read_cifti(SpaNov1_centrality_p2_file, brainstructures = "all")

# Correct again only across the two maps
new_correction <- correct_4_multiple_contrast_xiftis(list(SpaNov1_centrality_p1_xii, SpaNov1_centrality_p2_xii))

# Write to new cifti files
new_name <- str_replace(SpaNov1_centrality_p1_file, pattern = ".dscalar.nii", replacement = "new.dscalar.nii")
write_cifti(new_correction[[1]], cifti_fname = new_name)
## Writing left cortex.
## Writing right cortex.
## Writing subcortical data and labels.
## Creating CIFTI file from separated components.
new_name <- str_replace(SpaNov1_centrality_p2_file, pattern = ".dscalar.nii", replacement = "new.dscalar.nii")
write_cifti(new_correction[[2]], cifti_fname = new_name)
## Writing left cortex.
## Writing right cortex.
## Writing subcortical data and labels.
## Creating CIFTI file from separated components.
# Create data frame to compare differences
values <- c(as.matrix(SpaNov1_centrality_p1_xii), as.matrix(SpaNov1_centrality_p2_xii),
            as.matrix(new_correction[[1]]), as.matrix(new_correction[[2]]))
compare_p_dist <- data.frame(dist_type = rep(c("Correcting across main effect", "Only correcting across\npositive and negative"), each = length(values)/2),
                             log_p_values = values)

# Convert from log to actual p-values -log(0.05, 10)
compare_p_dist$p_values <- 10^-(compare_p_dist$log_p_values)

# Get lowest p-value
min_p <- min(compare_p_dist$p_values[compare_p_dist$dist_type == "Only correcting across\npositive and negative"])

# Five lowest p-values
low_5_p <- head(sort(compare_p_dist$p_values[compare_p_dist$dist_type == "Only correcting across\npositive and negative"]))

The five lowest p-values are 0.025789, 0.6326812, 0.8176687, 1, 1, 1.

Hippocampal gradient based on maximum statistic

# Create supplemental figure for maximum stat
## Creating a new data frame for combined plotting
HC_data_plotting2 <- HC_data_plotting
HC_data_agg_pos2$Hemisphere2 <- "Average"
HC_supp_fig_data <- data.frame(position = c(HC_data_plotting2$position, HC_data_agg_pos2$position),
                               max = c(HC_data_plotting2$max, HC_data_agg_pos2$max),
                               Hemisphere2 = c(HC_data_plotting2$Hemisphere2, HC_data_agg_pos2$Hemisphere2))

## Re-order Hemisphere2 so the order is left, right, then average
HC_supp_fig_data$Hemisphere2 <- factor(HC_supp_fig_data$Hemisphere2, 
                                       levels =  c("Left hippocampus", "Right hippocampus",
                                                   "Average"),
                                       ordered = TRUE)
## Create plot
scatter_max <- ggplot(HC_supp_fig_data, aes(x = -position, y = max)) + 
  facet_grid(~Hemisphere2, scales = "free_x") + 
  geom_point(aes(colour = max)) +
  geom_smooth(method = "lm", formula = y ~ x, colour = "black") +
  labs(x = "Distance to most anterior part (mm)", y = "Novelty Preference") + 
  scale_x_continuous(breaks =  seq(from = -45, to = 0, by = 15), labels = c("45", "30", "15", "0")) +
  scale_y_continuous(breaks = 1:6, labels = paste("Lvl", 1:6), limits = c(0.5, 6.5)) +
  scale_colour_viridis_c(option = "H", limits = c(1, 6), direction = -1) +
  theme_journal() +
  theme(strip.background = element_rect(color="white", fill="white")) +
  theme(legend.position = "none")

## Save with as 80 mm x 80 mm
ggsave(scatter_max, filename = paste0(figurePath, "Fig_S06_max_scatter.png"),
       dpi = dpi,  width = 120, height = 40, units = "mm")

Hippocampal gradient without averaging

# Create new copy of HC_data and create check sum
HC_data_no_avg <- HC_data
runCodeAgain3  <- check_if_md5_hash_changed(HC_data_no_avg, hash_table_name = "SpaNov_md5_hash_table.csv")

# Run/load the code/data
if(runCodeAgain3){
  # Seed
  set.seed(21312)
  
  # Other input
  lm_formula    <- "val ~ position * Hemisphere + tSNR + GS"
  nIter         <- 100000
  colName       <- "min"
  imageName     <- "SpaNov_permut_HC_grad_analysis1_cue-delay_no-avg.RData"
  
  grad_min_permut1_no_avg <- permutation_analysis(HC_data, lm_formula, nIter, 
                                                  colName, imageName)
  # Other input
  lm_formula    <- "val ~ position + tSNR + GS"
  nIter         <- 100000
  colName       <- "min"
  imageName     <- "SpaNov_permut_HC_grad_analysis1_cue-delay_L_no-avg.RData"
  
  grad_min_permut1_L_no_avg <- permutation_analysis(HC_data[HC_data$Hemisphere == "left", ], 
                                             lm_formula, nIter, colName, imageName)
  
  # Other input
  lm_formula    <- "val ~ position + tSNR + GS"
  nIter         <- 100000
  colName       <- "min"
  imageName     <- "SpaNov_permut_HC_grad_analysis1_cue-delay_R_no-avg.RData"
  
  grad_min_permut1_R_no_avg <- permutation_analysis(HC_data[HC_data$Hemisphere == 'right', ], 
                                             lm_formula, nIter, colName, imageName)
} else {
  load("intermediate_data/SpaNov_permut_HC_grad_analysis1_cue-delay_no-avg.RData")
  grad_min_permut1_no_avg   <- results
  load("intermediate_data/SpaNov_permut_HC_grad_analysis1_cue-delay_L_no-avg.RData")
  grad_min_permut1_L_no_avg <- results
  load("intermediate_data/SpaNov_permut_HC_grad_analysis1_cue-delay_R_no-avg.RData")
  grad_min_permut1_R_no_avg <- results
}

# Calculate p-value
dist      <- grad_min_permut1_no_avg$permuted_values[,5]
critVal   <- grad_min_permut1_no_avg$lm$coefficients[6]
grad_min_permut1_no_avg_p5 <- pValue_from_nullDist(critVal, dist, "two.sided")

# Calculate p-value
dist     <- grad_min_permut1_L_no_avg$permuted_values[,1]
critVal  <- grad_min_permut1_L_no_avg$lm$coefficients[2]
grad_min_permut1_L_no_avg_p1 <- pValue_from_nullDist(critVal, dist, "two.sided")

# Calculate p-value
dist      <- grad_min_permut1_R_no_avg$permuted_values[,1]
critVal   <- grad_min_permut1_R_no_avg$lm$coefficients[2]
grad_min_permut1_R_no_avg_p1 <- pValue_from_nullDist(critVal, dist, "two.sided")

# Calculate effect sizes
grad_min_permut1_no_avg_eff     <- eta_squared(car::Anova(grad_min_permut1_no_avg$lm, type = 2))
grad_min_permut1_L_no_avg_eff   <- eta_squared(car::Anova(grad_min_permut1_L_no_avg$lm, type = 2))
grad_min_permut1_R_no_avg_eff   <- eta_squared(car::Anova(grad_min_permut1_R_no_avg$lm, type = 2))

# Create data frame
## Create variables
tmp_formula   <- c("min ~ position * Hemisphere + tSNR + GS", "min (left) ~ position + tSNR + GS",
                   "min (right) ~ position + tSNR + GS")
tmp_coef_name <- c("interaction", "position", "position")
tmp_coef_val  <- c(grad_min_permut1_no_avg$lm$coefficients[6], grad_min_permut1_L_no_avg$lm$coefficients[2],
                   grad_min_permut1_R_no_avg$lm$coefficients[2])
tmp_coef_es   <- c(grad_min_permut1_no_avg_eff$Eta2_partial[5], grad_min_permut1_L_no_avg_eff$Eta2_partial[1],
                   grad_min_permut1_R_no_avg_eff$Eta2_partial[1])
tmp_coef_p    <- c(grad_min_permut1_no_avg_p5, grad_min_permut1_L_no_avg_p1, 
                   grad_min_permut1_R_no_avg_p1)

## Convert to numeric values and flip
tmp_coef_val <- -as.numeric(tmp_coef_val)

## Add to one data frame
tab_df <- data.frame(Formula = tmp_formula, Coefficient = tmp_coef_name,
                     beta = tmp_coef_val, eta_squared = tmp_coef_es, p = tmp_coef_p)

## Round columns
tab_df$beta        <- round(tab_df$beta, 3)
tab_df$eta_squared <- round(tab_df$eta_squared, 3)

## Create p-values by looping over all values
for(i in seq_along(tab_df$p)){
  tab_df$p[i] <-paste("p", pValue(as.numeric(tab_df$p[i])))
}

# Show table
kable(tab_df)
Formula Coefficient beta eta_squared p
min ~ position * Hemisphere + tSNR + GS interaction 0.016 0.004 p = .013
min (left) ~ position + tSNR + GS position 0.007 0.003 p = .158
min (right) ~ position + tSNR + GS position 0.026 0.032 p < .001
# Add better names for the hemispheres
HC_data$Hemisphere2 <- ifelse(HC_data$Hemisphere == "right",
                              "Right hippocampus", "Left hippocampus")

# Create plots
scatter_min_no_avg <- ggplot(HC_data, aes(x = -position, y = min)) + 
  facet_grid(~Hemisphere2, scales = "free_x") + 
  geom_jitter(aes(colour = min), height = 0.1, width = 0) +
  geom_smooth(method = "lm", formula = y ~ x, colour = "black") +
  labs(x = "Distance to most anterior part (mm)", y = "Novelty Preference") + 
  scale_x_continuous(breaks =  seq(from = -45, to = 0, by = 15), labels = c("45", "30", "15", "0")) +
  scale_y_continuous(breaks = 1:6, labels = paste("Lvl", 1:6), limits = c(0.5, 6.5)) +
  scale_colour_viridis_c(option = "H", limits = c(1, 6), direction = -1) +
  theme_journal() +
  theme(strip.background = element_rect(color="white", fill="white")) +
  theme(legend.position = "none")

## Save with as 80 mm x 80 mm
ggsave(scatter_min_no_avg, filename = paste0(figurePath, "Fig_S07_scatter_no-avg.png"),
       dpi = dpi,  width = 80, height = 40, units = "mm")

Relationship between spatial memory and successful memory encoding

# Calculate the average placement error for each object & subject. 
obj_avg <- ddply(OLM_7T_retrieval, c("subject", "objectName"), summarise,
                 euclideanDistance = mean(euclideanDistance), N = length(subject),
                 navTime = mean(navTime))

# Rank the objects within each subject
obj_avg <- ddply(obj_avg, c("subject"), mutate, rank = order(euclideanDistance))
obj_avg$f_rank <- factor(obj_avg$rank)

# Visualise results
p1 <- ggplot(obj_avg, aes(x = objectName, y = rank, fill = objectName)) + 
  geom_jitter(width = 0.25, height = 0, pch = 21) +
  geom_boxplot(outlier.shape = NA, alpha = 0.8) +
  stat_summary(geom = "point", fun = "mean", col = 'white', size = 0.8, shape = 24, aes(fill = objectName),
               position = position_dodge(width =  0.75),
               key_glyph = "rect") +
  scale_fill_viridis_d(option = "turbo") +
  labs(title = "Rank by objects",
       x = "Objects", y = "Within-subject rank")  +
  theme_journal() +
  theme(legend.position = 'none',
        plot.title = element_text(hjust = 0.5))

p2 <- ggplot(obj_avg, aes(x = f_rank, y = euclideanDistance, fill = f_rank)) + 
  geom_jitter(width = 0.25, height = 0, pch = 21) +
  geom_boxplot(outlier.shape = NA, alpha = 0.8) +
  stat_summary(geom = "point", fun = "mean", col = 'white', size = 0.8, shape = 24, aes(fill = f_rank),
               position = position_dodge(width =  0.75),
               key_glyph = "rect") +
  labs(title = "Placement error by rank",
       x = "Rank", y = "Avg. placement error (vm)")  +
  theme_journal() +
  scale_fill_viridis_d() +
  theme(legend.position = 'none',
        plot.title = element_text(hjust = 0.5))

p3 <- ggplot(obj_avg, aes(x = f_rank, y = navTime, fill = f_rank)) + 
  geom_jitter(width = 0.25, height = 0, pch = 21) +
  geom_boxplot(outlier.shape = NA, alpha = 0.8) +
  stat_summary(geom = "point", fun = "mean", col = 'white', size = 0.8, shape = 24, aes(fill = f_rank),
               position = position_dodge(width =  0.75),
               key_glyph = "rect") +
  labs(title = "Navigation time by rank",
       x = "Rank", y = "Navigation time (seconds)")  +
  theme_journal() +
  scale_fill_viridis_d() +
  theme(legend.position = 'none',
        plot.title = element_text(hjust = 0.5))

###############################################################################_
# Create base path
base_path <- "/media/alex/work/Seafile/imaging_results/SpaNov/OLMe_7T_encoding1_cue-delay_smo4_MSMAll/cope7.feat/stats/vwc/"

# Load the beta values for encoding contrast
encoding_betaMap_xii  <- read_cifti(paste0(base_path, "Y1.dtseries.nii"), brainstructures = "all")

# Get p-value maps
encoding_FDRp_c1_xii <- read_cifti(paste0(base_path, "results_lvl2cope1_dat_ztstat_cfdrp_c1.dscalar.nii"), brainstructures = "all")
encoding_FDRp_c2_xii <- read_cifti(paste0(base_path, "results_lvl2cope1_dat_ztstat_cfdrp_c2.dscalar.nii"), brainstructures = "all")
# Count significant gray-ordinates
# sum(as.matrix(encoding_FDRp_c1_xii) > cutOff) # = 548
# sum(as.matrix(encoding_FDRp_c2_xii) > cutOff) # = 548

# Calculate average successful encoding contrast in novelty masls
encoding_WB_familiarity <- as.matrix(encoding_betaMap_xii)[as.matrix(WB_familiarity_mask) == 1, ]
encoding_WB_familiarity <- colMeans(encoding_WB_familiarity)
encoding_WB_novelty     <- as.matrix(encoding_betaMap_xii)[as.matrix(WB_novelty_mask) == 1, ]
encoding_WB_novelty     <- colMeans(encoding_WB_novelty)

# Calculate stats
x      <- encoding_WB_familiarity
d1     <- signif(mean(x)/sd(x), digits1)
str1   <- mean_SD_str2(x, report_type, digits1, rounding_type)
BF10_1 <- reportBF(ttestBF(x))
x      <- encoding_WB_novelty
d2     <- signif(mean(x)/sd(x), digits1)
str2   <- mean_SD_str2(x, report_type, digits1, rounding_type)
BF10_2 <- reportBF(ttestBF(x))

# Create a data frame for plotting
encoding_WB_df <- data.frame(mask = rep(c("Familiarity", "Novelty"), each = length(encoding_WB_familiarity)),
                             contrast = c(encoding_WB_familiarity, encoding_WB_novelty))

p4 <- ggplot(encoding_WB_df, aes(x = mask, y = contrast, fill = mask)) + 
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_jitter(width = 0.25, height = 0, pch = 21) +
  geom_boxplot(outlier.shape = NA, alpha = 0.8) +
  stat_summary(geom = "point", fun = "mean", col = 'white', size = 0.8, shape = 24, aes(fill = mask),
               position = position_dodge(width =  0.75),
               key_glyph = "rect") +
  annotate("text", x = 1, y = 200, label = paste("BF10 =", BF10_1), size = 2) +
  annotate("text", x = 2, y = 250, label = paste("BF10 =", BF10_2), size = 2) +
  labs(title = "",
       x = "Mask", y = "Encoding success (a.u.)")  +
  theme_journal() + 
  scale_fill_manual(values = cool_warm_colours) +
  theme(legend.position = 'none',
        plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5),
        plot.margin = unit(c(0, 1, 0, 22.5), "mm"))

# Assemble the plot
p <- plot_grid(p1, p2, p3, NULL, NULL, p4,
               align = "h", nrow = 2, #labels = c("a", "b", "c", "d", "", "e"), 
               label_size = 10, rel_heights = c(1, 1))
# Save
ggsave(p,
       filename = paste0(figurePath, "Fig_S08_successful_encoding.png"), 
       dpi = dpi,
       width = 180,
       height = 100,
       units = "mm")
  • Encoding effect in familiarity regions: d = -0.068, M = -2.9 (SD = 43), BF10 = 0.16
  • Encoding effect in novelty regions: d = 0.43, M = 31 (SD = 73), BF10 = 12.99

Distribution of principal functional connectivity gradient scores across spatial novelty preference levels in the posterior medial cortex

# Compare the PMC gradient with the first connectivity gradient from Margulies et al. (2016)
## Load PMC gradient
PMC_gradient <- read_cifti("cifti_results/SpaNov_gradient_PMC_min_cue-delay.dlabel.nii")

## Add to Margulies_grad data frame
Margulies_grad$PMC_gradient <- c(PMC_gradient$data$cortex_left, PMC_gradient$data$cortex_right)

## Subset to only vertices inside the PMC gradient
Margulies_grad_PMC <- Margulies_grad[Margulies_grad$PMC_gradient != 0, ]

## Make factor
Margulies_grad_PMC$PMC_gradient <- factor(Margulies_grad_PMC$PMC_gradient, levels = 1:6, ordered = TRUE)

## Remove Level5 because there are only 3 vertices
Margulies_grad_PMC <- Margulies_grad_PMC[Margulies_grad_PMC$PMC_gradient != 5, ]

## Create a figure
Margulies_novelty <- ggplot(Margulies_grad_PMC, aes(x = PMC_gradient, y = Grad_1, fill = PMC_gradient)) +
  geom_jitter(height = 0, size = 0.5, alpha = 0.1, width = 0.2) +
  geom_boxplot(alpha = 0.8, outlier.shape = NA) +
  stat_summary(geom = "point", fun = "mean", col = 'white', shape = 24, 
               position = position_dodge(width =  0.75), size = 0.8) +
  scale_fill_manual(values = novFam_gradient[-5]) +
  theme_journal() + 
  theme(legend.position = "none") +
  labs(x = "Spatial novelty-familarity gradient (PMC)", y = TeX(r"( $Margulies'\ et\ al^{10}\ 1^{st}\ gradient$ )"))

# Save
ggsave(Margulies_novelty,
       filename = paste0(figurePath, "Fig_S09_PMC_gradient_Margulies.png"), 
       dpi = dpi,
       width = 80,
       height = 60,
       units = "mm")

Hippocampal connectivity patterns correlating with Margulies gradient

# Load results
load("data/ignore_fMRI_version1/grad_FRSC/HC_2_cortex_FRSC_Margulies_gradient.RData")

# Apply fisher z-transformation
HC_2_cortex_Margulies_gradient$z_correlation <-  z_transform_fisher(HC_2_cortex_Margulies_gradient$correlation)

# Add MNI y coordinate to it 
HC_2_cortex_Margulies_gradient$MNI_y <- rep(rep(HC_data$y, 10), 56)

# Divide into anterior and poster
HC_2_cortex_Margulies_gradient$AP   <- ifelse(HC_2_cortex_Margulies_gradient$MNI_y >= -21, 
                                              "Anterior", "Posterior")
# Calculate average
HC_2_Margulies_agg1 <- ddply(HC_2_cortex_Margulies_gradient, c("Rnum", "Gradient", "AP"),
                            summarise, z_correlation = mean(z_correlation))

# Make gradient a factor
HC_2_Margulies_agg1$Gradient <- factor(HC_2_Margulies_agg1$Gradient)

# Calculate differences
HC_2_Margulies_agg2 <- ddply(HC_2_Margulies_agg1, c("Rnum", "Gradient"),
                            summarise, diff = mean(z_correlation[1] - z_correlation[2]))

# Calculate tests for each gradient
BF10 <- rep(NA, 10)
d    <- rep(NA, 10)

for(i in 1:10){
  # Calculate one sample difference
  diff <- HC_2_Margulies_agg2$diff[HC_2_Margulies_agg2$Gradient == i]
  d[i] <- mean(diff)/sd(diff)
  BF10[i] <- signif(reportBF(ttestBF(diff)), 2)
}

# Make it a data frame for simpler adding
annot_df <- data.frame(label = paste("BF10 =", BF10),
                       Gradient = factor(1:10),
                       y = rep(c(0.09, 0.08), 5))

# Visualise
HC_2_margulies <- ggplot(HC_2_Margulies_agg2, aes(x = Gradient, y = diff, fill = Gradient)) +
  geom_hline(yintercept = 0, linetype = "dashed", colour = "darkgrey") +
  scale_fill_viridis_d() +
  geom_boxplot(outlier.shape = NA) + 
  stat_summary(geom = "point", fun = "mean", col = 'white', shape = 24, 
               position = position_dodge(width =  0.75), size = 0.8) +
  annotate("text", x = annot_df$Gradient, y = annot_df$y, label = annot_df$label, size = 1.5) +
  theme_journal() +
  theme(legend.position = "none", plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5)) +
  labs(x = "Principal gradient", y = "Anterior - posterior difference\nin correlation")

# Save
ggsave(HC_2_margulies,
       filename = paste0(figurePath, "Fig_S10_HC_2_Margulies.png"), 
       dpi = dpi,
       width = 120,
       height = 60,
       units = "mm")

Posterior predictive checks for the novelty score model

Find information here: https://mc-stan.org/rstanarm/reference/pp_check.stanreg.html

# Set extra seed because these plot are relatively variable with regard to the tails
set.seed(20240206)

p1 <- brms::pp_check(m_noveltyScore_run1, ndraws = 100) + 
  scale_colour_manual(values = c("black", base_col[1]), name = "") +
  coord_cartesian(xlim = c(-20, 20), expand = TRUE) + 
  theme_journal() + 
  theme(axis.line.y = element_blank(),
        axis.text.y = element_blank(), 
        axis.ticks.y = element_blank(),
        legend.position = "none") 

p2 <- brms::pp_check(m_noveltyScore_run2, ndraws = 100) + 
  scale_colour_manual(values = c("black", base_col[1]), name = "") +
  coord_cartesian(xlim = c(-50, 50), expand = TRUE) + 
  theme_journal() + 
  theme(axis.line.y = element_blank(),
        axis.text.y = element_blank(), 
        axis.ticks.y = element_blank(),
        legend.position = "none") 

p3 <- brms::pp_check(m_conn2, ndraws = 100) + 
  scale_colour_manual(values = c("black", base_col[1]), name = "") +
  coord_cartesian(xlim = c(-10, 10), expand = TRUE) + 
  theme_journal() + 
  theme(axis.line.y = element_blank(),
        axis.text.y = element_blank(), 
        axis.ticks.y = element_blank(),
        legend.position = "none") 

p4 <- brms::pp_check(m_navTime3, ndraws = 100) + 
  scale_colour_manual(values = c("black", base_col[1]), name = "") +
  coord_cartesian(xlim = c(0, 60), expand = TRUE) + 
  theme_journal() + 
  theme(axis.line.y = element_blank(),
        axis.text.y = element_blank(), 
        axis.ticks.y = element_blank(),
        legend.position = "none") 

p6 <- brms::pp_check(m_visits_run1, ndraws = 100) + 
  scale_colour_manual(values = c("black", base_col[1]), name = "") +
  coord_cartesian(xlim = c(0, 20), expand = TRUE) + 
  theme_journal() + 
  theme(axis.line.y = element_blank(),
        axis.text.y = element_blank(), 
        axis.ticks.y = element_blank(),
        legend.position = "none") 

p7 <- brms::pp_check(m_visits_run2, ndraws = 100) + 
  scale_colour_manual(values = c("black", base_col[1]), name = "") +
  coord_cartesian(xlim = c(0, 20), expand = TRUE) + 
  theme_journal() + 
  theme(axis.line.y = element_blank(),
        axis.text.y = element_blank(), 
        axis.ticks.y = element_blank(),
        legend.position = "none") 

p8 <- brms::pp_check(m_timeSince_run1, ndraws = 100) + 
  scale_colour_manual(values = c("black", base_col[1]), name = "") +
  coord_cartesian(xlim = c(0, 500), expand = TRUE) + 
  theme_journal() + 
  theme(axis.line.y = element_blank(),
        axis.text.y = element_blank(), 
        axis.ticks.y = element_blank(),
        legend.position = "none") 

p9 <- brms::pp_check(m_timeSince_run2, ndraws = 100) + 
  scale_colour_manual(values = c("black", base_col[1]), name = "") +
  coord_cartesian(xlim = c(0, 2000), expand = TRUE) + 
  theme_journal() + 
  theme(axis.line.y = element_blank(),
        axis.text.y = element_blank(), 
        axis.ticks.y = element_blank(),
        legend.position = "none") 

# Solution to change the linewidth found here: https://stackoverflow.com/questions/55450441/how-to-change-size-from-specific-geom-in-ggplot2
p1$layers[[1]]$aes_params$linewidth <- 0.1
p1$layers[[2]]$aes_params$linewidth <- 0.1
p2$layers[[1]]$aes_params$linewidth <- 0.1
p2$layers[[2]]$aes_params$linewidth <- 0.1
p3$layers[[1]]$aes_params$linewidth <- 0.1
p3$layers[[2]]$aes_params$linewidth <- 0.1
p4$layers[[1]]$aes_params$linewidth <- 0.1
p4$layers[[2]]$aes_params$linewidth <- 0.1
p6$layers[[1]]$aes_params$linewidth <- 0.1
p6$layers[[2]]$aes_params$linewidth <- 0.1
p7$layers[[1]]$aes_params$linewidth <- 0.1
p7$layers[[2]]$aes_params$linewidth <- 0.1
p8$layers[[1]]$aes_params$linewidth <- 0.1
p8$layers[[2]]$aes_params$linewidth <- 0.1
p9$layers[[1]]$aes_params$linewidth <- 0.1
p9$layers[[2]]$aes_params$linewidth <- 0.1

# Combine plots
p_comb <- plot_grid(p1, p2, p3, p4, p6, p7, p8, p9, 
                    labels = "auto", label_size = 10, align = "hv")

ggsave(p_comb,
       filename = paste0(figurePath, "Fig_S11_ppc.png"), 
       dpi = dpi,
       width = 90,
       height = 90,
       units = "mm")

Distribuion of starting relative to object locations

# Use one participant for visualisation because they are all the same anyway
start_vis_data <- OLM_7T_encoding[OLM_7T_encoding$subject == unique(OLM_7T_encoding$subject)[1], ]

# Remove the gift box that location was variable
no_gift <- obj_locations[obj_locations$objectName != "gift", ]

# Create the plot
start_locations <- ggplot() +
  # Circles
  geom_polygon(data = circle1, mapping = aes(x = x, y = z), fill = "lightgrey") +
  geom_polygon(data = circle2, mapping = aes(x = x, y = z), fill = "darkgrey") +
  geom_point(data = start_vis_data, mapping = aes(x = start_x, y = start_z, colour = targetNames)) +
  geom_segment(data = start_vis_data, mapping = aes(x = start_x, y = start_z, xend = object_x, yend = object_z, 
                                                    colour = targetNames)) +
  geom_point(data = no_gift, mapping = aes(x = object_x, y = object_z, colour = targetNames), size = 2) +
  scale_colour_viridis_d(option = "turbo") +
  # Other settings
  coord_equal(xlim = c(-90, 90)) +
  theme_void() +
  theme(plot.background = element_rect(fill = "white", colour = "white")) +
  labs(colour = "") +
  theme(legend.text = element_text(size=rel(0.5)),
        legend.key.height = unit(0.5,"line"),
        legend.key.width = unit(0.1,"line"),
        legend.title = element_blank(),
        legend.box.margin = margin(0, 0, 0, -10))

ggsave(start_locations,
       filename = paste0(figurePath, "Fig_S12_start_locations.png"), 
       dpi = dpi, width = 80, height = 80, units = "mm")

Correlation between novelty measures

###############################################################################_
# Number of visits
# Set range
x1_range   <- range(m_visits_run1$data$s_onset_rel)
x1_points  <- seq_min_2_max(m_visits_run1$data$s_onset_rel)
x2_range   <- range(m_visits_run2$data$s_onset_rel) 
x2_points  <- seq_min_2_max(m_visits_run2$data$s_onset_rel)

# Get subject-level predictions
pred1_df <- predictions(m_visits_run1, newdata = datagrid(s_onset_rel = x1_points, subject = unique), 
                           by = c("s_onset_rel", "subject"))
pred2_df <- predictions(m_visits_run2, newdata = datagrid(s_onset_rel = x2_points, subject = unique), 
                           by = c("s_onset_rel", "subject"))
pred1_df <- as.data.frame(pred1_df)
pred2_df <- as.data.frame(pred2_df)

# Calculate the minimum for each subject, which is used for colouring
pred1_df <- ddply(pred1_df, c("subject"), mutate, min = min(estimate))
pred2_df <- ddply(pred2_df, c("subject"), mutate, min = min(estimate))

# Create plots
## Run 1
p1_1 <- ggplot(data = pred1_df, aes(x = s_onset_rel, y = estimate, group = subject, colour = min)) +
  geom_line() +
  scale_color_viridis_c() + theme_journal() +
  theme(legend.position = "none",
        plot.margin = unit(c(1, 2, 1, 2), "mm")) +
  labs(title = "Run 1", x = "Time", y = "Number of visits") +
  scale_x_continuous(breaks = x1_range, labels = c("Start", "End")) +
  scale_y_continuous(breaks = c(2, 6, 10)) +
  coord_cartesian(xlim = x1_range, ylim = c(2, 10), expand = FALSE) 

## Run 2
p1_2 <- ggplot(data = pred2_df, aes(x = s_onset_rel, y = estimate, group = subject, colour = min)) +
  geom_line() +
  scale_color_viridis_c() + theme_journal() +
  theme(legend.position = "none",
        plot.margin = unit(c(1, 2, 1, 2), "mm")) +
  labs(title = "Run 2", x = "Time", y = "Number of visits") +
  scale_x_continuous(breaks = x2_range, labels = c("Start", "End")) +
  scale_y_continuous(breaks = c(6, 10, 14)) +
  coord_cartesian(xlim = x2_range, ylim = c(6, 14), expand = FALSE) 

###############################################################################_
# Get original scale for time since
# Load 
load("event_tables/images/SpaNov_event_file_gradients.RData")

# Convert list to data frame
data <- rbindlist(subj_list, idcol = "subject")

# Add subjects' R number
data$ppid <- subjIDs_R[data$subject]

# Function to match runStartTime
find_runStartTime <- function(ppid, run){
  # Get corresponding to find the run in the trial data
  anonKey <- lookupTable$anonKey[lookupTable$Rnum == ppid[1]]
  
  # Use the anonKey & run to get runStartTime
  runStartTime <- OLM_7T_trial_results$runStartTime[OLM_7T_trial_results$subject == anonKey & 
                                                    OLM_7T_trial_results$block_num == run[1]]
  
  return(runStartTime[1])
}

data <- ddply(data, c("subject", "ppid", "run"), mutate, 
              runStartTime = find_runStartTime(ppid, run))

data$subject <- as.character(data$subject)

# Make time relative to the start of the run. The real onset times will be 
# slightly different but this will not matter for this. 
data$onset_rel <- data$onset - data$runStartTime

# Add run type
data$runType <- "encoding"
data$runType[data$run == 2 | data$run == 4] <- "retrieval"

# Subset to only encoding
data_sub <- data[data$runType == "encoding", ]

# Change run number to match the description in paper. Run 3 is Encoding Run 2
data_sub$run[data_sub$run == 3] <- 2

# Convert run to factor
data_sub$f_run <- as.factor(data_sub$run)

# Run 1 
## Remove NA and get correct run
## Use unique name for data frame because of marginaleffects issues
m_timeSince_run1_data <- na.omit(data_sub[data_sub$run == 1, ])

## Scale x and y
m_timeSince_run1_data$s_onset_rel  <- scale_m0_sd0.5(m_timeSince_run1_data$onset_rel)
m_timeSince_run1_data$s_timeSince  <- drop(scale(m_timeSince_run1_data$lastVisit)) 
# Use drop() to deal with marginal effects issues

# Run 2 
## Remove NA and get correct run
## Use unique name for data frame because of marginaleffects issues
m_timeSince_run2_data <- na.omit(data_sub[data_sub$run == 2, ])

## Scale x and y
m_timeSince_run2_data$s_onset_rel  <- scale_m0_sd0.5(m_timeSince_run2_data$onset_rel)
m_timeSince_run2_data$s_timeSince  <- drop(scale(m_timeSince_run2_data$lastVisit))
# Use drop() to deal with marginal effects issues

###############################################################################_
# Time since
# Set range
x1_range   <- range(m_timeSince_run1$data$s_onset_rel)
x1_points  <- seq_min_2_max(m_timeSince_run1$data$s_onset_rel)
x2_range   <- range(m_timeSince_run2$data$s_onset_rel) 
x2_points  <- seq_min_2_max(m_timeSince_run2$data$s_onset_rel)

# Get subject-level predictions
pred1_df <- predictions(m_timeSince_run1, newdata = datagrid(s_onset_rel = x1_points, subject = unique), 
                           by = c("s_onset_rel", "subject"))
pred2_df <- predictions(m_timeSince_run2, newdata = datagrid(s_onset_rel = x2_points, subject = unique), 
                           by = c("s_onset_rel", "subject"))
pred1_df <- as.data.frame(pred1_df)
pred2_df <- as.data.frame(pred2_df)

# Calculate the minimum for each subject, which is used for colouring
pred1_df <- ddply(pred1_df, c("subject"), mutate, min = min(estimate))
pred2_df <- ddply(pred2_df, c("subject"), mutate, min = min(estimate))

# Create plots
## Run 1
p2_1 <- ggplot(data = pred1_df, aes(x = s_onset_rel, y = estimate, group = subject, colour = min)) +
  geom_line() +
  scale_color_viridis_c() + theme_journal() +
  theme(legend.position = "none",
        plot.margin = unit(c(1, 2, 1, 2), "mm")) +
  labs(title = "Run 1", x = "Time", y = "Time elapsed (seconds)") +
  scale_x_continuous(breaks = x1_range, labels = c("Start", "End")) +
  scale_y_continuous(breaks = c(50, 125, 200)) +
  coord_cartesian(xlim = x1_range, ylim = c(50, 200), expand = FALSE) 

## Run 2
p2_2 <- ggplot(data = pred2_df, aes(x = s_onset_rel, y = estimate, group = subject, colour = min)) +
  geom_line() +
  scale_color_viridis_c() + theme_journal() +
  theme(legend.position = "none",
        plot.margin = unit(c(1, 2, 1, 2), "mm")) +
  labs(title = "Run 2", x = "Time", y = "Time elapsed (seconds)") +
  scale_x_continuous(breaks = x2_range, labels = c("Start", "End")) +
  scale_y_continuous(breaks = c(50, 300, 550)) +
  coord_cartesian(xlim = x2_range, ylim = c(50, 550), expand = FALSE) 

###############################################################################_
# Calculate correlation
temp_corr <- cor.test(event_info_df_encoding$visits, event_info_df_encoding$lastVisit)
plot_str  <- paste("r", rValue(temp_corr$estimate), "p", pValue(temp_corr$p.value))

# Create scatter plot
corr_novelty_measures <- ggplot(event_info_df_encoding, aes(x = visits, y = lastVisit)) +
  geom_point(alpha = 0.1) + 
  geom_smooth(method = "lm", formula = "y ~ x", colour = base_col[1]) +
  annotate("text", x = I(0.8), y = I(0.8), label = plot_str, size = 2) +
  theme_journal() +
  theme(legend.position = "none",
        plot.title = element_text(hjust = 0.5)) +
  labs(x = "Number of visits", y = "Time elapsed (seconds)",
       title = "Relationship between novelty measures")

###############################################################################_
# Combine into one figure
combined_plot <- plot_grid(corr_novelty_measures, plot_grid(p1_1, p1_2, p2_1, p2_2, 
                                                            ncol = 2, labels = c("b", "", "c", ""), 
                                                            align = "hv",  label_size = 10),
                           ncol = 2, labels = c("a", ""), label_size = 10)

# Save
ggsave(combined_plot,
       filename = paste0(figurePath, "Fig_S13_corr_novelty_measures.png"), 
       dpi = dpi, width = 190, height = 82, units = "mm")